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A  new  three-dimensional  explicit  Navier-Stokes  procedure  has  been  developed  for 
computation  of  turbulent  turbomachinery  flows.  Several  numerical  strategies  and 
modelling  techniques  have  been  developed  and  incorporated  which  enable  convergent  and 
accurate  predictions  of  high  Reynolds  number  flowfields  across  a  wide  range  of  Mach 
numbers.  These  include  incorporation  of  a  compressible  low  Reynolds  number  form  of  the 
k-t  turbulence  model  in  a  fully  explicit  fashion,  appropriate  stability  bound  treatment  of  the 
transport  turbulence  model  and  other  physical  and  solution  parameters,  eigenvalue  and  local 
velocity  artificial  dissipation  scalings,  a  compact  flux  evaluation  procedure  and  a  hybrid 
low  Reynolds  number  k-e  /  algebraic  Reynolds  stress  model. 

The  derivation  of  the  governing  transport  equations  and  the  numerical 
discretization  procedures  are  provided  for  three  dimensions,  although  two-dimensional 
computations  are  also  considered.  A  compact  central  finite  difference  flux  evaluation 
scheme,  is  introduced.  This  scheme  allows  for  the  straightforward  application  of  boundary 
conditions  often  attributed  to  finite  volume  formulations,  and  eliminates  the  metric 
singularity  at  the  interface  between  solid  and  periodic  boundaries. 

Detailed  stability  and  order  of  magnitude  analyses  are  performed  on  the  discrete 
system  of  seven  governing  equations.  The  results  of  these  analyses  and  corroborative 
numerical  experiments  are  provided.  Conclusions  are  drawn  concerning  the  influence  of 
system  rotation  and  turbulence  transport  source  terms,  grid  clustering,  effective  diffusivity, 
artificial  dissipation,  convective  acceleration  terms,  implicit  source  term  treatment  and  the 
coupling  of  the  discrete  mean  flow  equation  system  to  the  turbulence  model  equations  on 
the  stability  of  the  numerical  scheme.  Notable  is  the  conclusion  that,  contrary  to  general 
perception,  turbulence  model  source  terms  themselves  do  not  adversely  affect  the  stability 


of  the  numerical  scheme  used,  for  many  low  Reynolds  number  forms,  provided  that  certain 
constraints  are  met.  Applicability  of  these  results  to  explicit  multigrid  and  implicit  time 
marching  schemes  is  discussed. 

The  highly  stretched  grids,  required  to  resolve  near  wall  flow  physics  in  high 
Reynolds  number  turbomachinery  flow  computations,  give  rise  to  important  artificial 
dissipation  issues.  Accordingly,  alternative  local  velocity  and  flux  Jacobian  eigenvalue 
scalings  are  provided  for  the  dissipation  operators.  It  is  illustrated  that  such  scalings  are 
required  for  convergent  and  accurate  solutions  of  these  flows. 

A  hybrid  k-e  /  algebraic  Reynolds  stress  model  (ARSM)  has  been  developed.  This 
hybrid  approach  reconciles  the  near  wall  damping  provided  by  the  low  Reynolds  number  k- 
e  model  with  the  high  Reynolds  number  form  ARSM.  The  accuracy  of  this  approach  is 
validated  for  a  simple  shear  layer,  and  the  numerical  stability  of  the  model  is  analyzed  and 
verified. 

Two-dimensional  results  are  provided  for  a  supersonic  compressor  cascade 
operating  at  unique  incidence  condition  and  a  low-subsonic  double  circular  arc  compressor 
cascade.  Detailed  comparison  of  the  numerical  solution  with  experimental  results  is 
provided.  Shock  structure,  boundary  layer  and  wake  physics  are  well  predicted,  providing 
validation  of  the  numerical  and  turbulence  modelling  approaches. 

Three-dimensional  validation  is  provided  by  the  results  of  an  incompressible 
curved  duct  flow  computation,  for  which  wall  static  pressure  and  LDA  primary  and 
secondary  velocity  measurements  are  available.  A  high  Reynolds  number  axial  rotor  flow, 
for  which  extensive  experimental  data  is  also  available,  was  computed.  This  large  scale 
computation  is  shown  to  capture  rotational  inviscid,  blade  and  endwall  boundary  layer  and 
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wake  flow  physics,  including  spanwise  mixing  effects,  with  good  accuracy. 

Finally,  a  backswept  transonic  centrifugal  compressor  flow,  for  which  L2F 
meridional  passage  velocity  measurements  are  available,  is  computed.  Such  flows  are  the 
most  complex  that  exist  in  turbomachines.  Full  Navier-Stokes  solutions  are  presented 
which  are  shown  to  capture  detailed  viscous  dominated  flow  features,  including  tip 
clearance  and  curvature  induced  and  rotation  induced  secondary  motions,  with  good 
accuracy.  Relative  helicity  is  used  to  help  interrogate  these  computational  results,  and  it  is 
found  that  this  parameter  is  very  useful  in  providing  insight  into  the  secondary  motions  in 
these  machines.  Results  of  impeller  flowfield  calculations  using  the  hybrid  model  with  and 
without  the  influence  of  system  rotation  are  also  provided  and  interpreted. 
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CHAPTER  1 

INTRODUCTION 


Computation  of  viscous  flows  by  numerical  solution  of  the  Navier-Stokes 
equations  has  become  increasingly  feasible  due  to  improvements  in  numerical  algorithms, 
and  to  the  ever-increasing  speed  and  memory  of  digital  computers.  The  most  up  to  date 
CFD  codes  available  today  are  capable  of  calculating  steady  3-D  viscous  flows  about  entire 
vehicles,  and  even  unsteady  viscous  flows  in  3-D  turbomachinery  stages.  However, 
despite  the  rapid  advance  towards  exploiting  the  power  of  computers  now  available,  some 
serious  limitations  of  these  codes  have  yet  to  be  adequately  resolved.  Surely  the  most 
profound  of  these  limitations  is  the  lack  of  accurate,  general  turbulence  models.  Secondary 
to  this,  but  of  much  concern,  is  the  role  of  artificial  dissipation  in  Navier-Stokes 
calculations. 

The  flowfield  in  all  turbomachines  is  characterized  by  significant  rotational  inviscid 
phenomena,  arising  from  stagnation  enthalpy  and  pressure  gradients  due  to  shaft  work, 
flow  turning  and  compressibility  effects.  Consequently,  a  great  deal  of  the  physics  in  these 
machines  can  be  analyzed  with  inviscid  methods.  For  this  reason,  inviscid  solution 
procedures  have  advanced  to  the  point,  and  are  efficient  enough,  that  they  are  today  used 
routinely  in  turbomachinery  design. 

Alternatively,  viscous  flow  phenomena  play  a  dominant  role  in  many  parameters  of 
interest  to  the  engineer.  Specifically,  machine  performance  parameters  such  as  efficiency, 
maximum  pressure  rise/drop  and  attainable  mass  flow  rates  are  strongly  influenced  by 
viscous  effects.  The  developing  boundary  layers  on  blade  and  endwall  surfaces  give  rise  to 
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losses  and  blockage.  The  convection  and  diffusion  of  blade  row  wakes  and  some 
spanwise  mixing  effects  are  viscous  phenomena  which  play  important  roles  in  multistage 
machines.  Also  of  large  concern,  especially  in  turbines  and  centrifugal  machines,  are 
secondary  motions  which  arise  due  to  the  turning  and  rotation  of  shear  layers  in  these 
machines.  The  direct  effect  of  apparent  stresses,  due  to  flowfield  turbulence,  is  a  factor  in 
all  turbomachines  since  Reynolds  numbers  are  characteristically  large.  The  physics  of  tip 
clearance  vortices,  localized  flow  reversal,  rotating  stall,  surge  and  choking  are  all 
intimately  tied  to  viscous  flow  effects.  In  order  to  develop  a  further  understanding  of  these 
phenomena,  nothing  short  of  analyses  that  incorporate  the  full  Reynolds  averaged  Navier- 
Stokes  equations  are  sufficient 

The  engineer  needs  analysis  tools  which  are  1)  accurate,  2)  versatile,  3)  relatively 
easy  to  use  on  a  routine  basis,  and  4)  fast/efficient  Unfortunately,  for  these  very  reasons, 
the  current  state  of  the  art  does  not  allow  for  detailed  Navier-Stokes  analysis  in  the 
turbomachinery  design  environment.  At  present,  high  Reynolds  number,  three- 
dimensional  turbomachinery  passage  flows  have  been  computed  by  a  number  of 
researchers  using  on  the  order  of  105  grid  points,  but  the  central  memory  and  CPU  run 
time  requirements  for  such  calculations  are  large.  Also,  even  grids  of  this  size  are  not 
adequate  to  provide  grid  converged,  detailed  viscous  flow  predictions.  The  state  of  the  art 
seems  to  be,  rather,  the  reasonable  prediction  of  global  three-dimensional  vorticity 
distributions  at  design  and  off  design  conditions.  Also,  in  addition  to  accuracy  limitations 
imposed  by  practical  grid  size  constraints,  simple  turbulence  models  and  the  inescapable 
presence  of  artificial  dissipation  introduce  additional  inaccuracies. 

Pre-processing,  including  grid  generation  and  flow  initialization,  can  be  daunting 
and  time-consuming  tasks  for  complex  turbomachinery  configurations.  Application  of 
Navier-Stokes  flow  solvers  themselves  generally  require  an  experienced  analyst,  with 
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detailed  understanding  of  the  physical  models  and/or  numerical  procedures  incorporated  in 
these  codes.  Lastly,  postprocessing  large  scale  solutions  can  take  days  of  an  analyst’s  time 
to  extract  useful  design  information.  So,  in  short,  Navier-Stokes  methodologies  currently 
fail  on  several  counts,  from  an  engineer's  standpoint.  For  this  reason,  there  is  currently  a 
great  deal  of  research  to  improve  the  usefulness  of  these  methods,  with  emphasis  on 
turbulence  and  other  physical  modeling,  incorporation  of  efficient  algorithms  and  improved 
computer  architectures,  and  the  use  of  full  scale  simulation  studies  to  develop  an 
understanding  of  particular  physical  processes  in  turbomachines  [see  recent  paper  by 
Leylek  and  Wisler  (1990),  for  an  good  example  of  this  last  item)]. 

The  purpose  of  this  thesis  is  to  develop  a  new  explicit,  three-dimensional  Navier- 
Stokes  methodology  for  turbomachinery  flowfield  prediction.  In  accordance  with  the  need 
to  address  the  items  mentioned  above,  emphasis  has  been  placed  on  1)  incorporation  of  a 
low  Reynolds  number  two-equation  turbulence  model  in  an  explicit  numerical  solution 
procedure,  2)  comprehensive  investigation  of  the  stability  of  the  procedure,  including  the 
influence  of  various  geometric  and  flowfield  parameters  on  numerical  stability,  with 
emphasis  on  the  turbulence  model,  3)  investigation  of  the  role  of  artificial  dissipation  in  the 
accuracy  and  stability  of  viscous  flow  computations  on  highly  stretched  grids,  4) 
simulation  of  a  variety  of  complex  two-  and  three-dimensional  turbomachinery  flows, 
across  a  wide  range  of  Mach  numbers,  including  a  transonic  centrifugal  compressor,  with 
the  goal  of  using  the  flow  solution  to  provide  a  detailed  understanding  of  the  viscous 
physical  mechanisms  in  these  machines. 

In  three-dimensional  Navier-Stokes  calculations  of  steady  state  turbomachinery 
flowfields,  noth  explicit  and  implicit  time-marching  procedures  have  been  implemented  by 
a  number  of  researchers.  As  mentioned  above,  high  Reynolds  number  turbomachinery 
passage  flows  have  been  computed  by  a  number  of  researchers  using  on  the  order  of  105 
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grid  points  [See  recent  publications  by  Bansod  and  Rhie  (1990),  Chima  and  Yokota 
(1988),  Choi  and  Knight  (1989),  Walker  and  Dawes  (1989),  Hah  and  Wennerstrom 
(1990),  Hobson  and  Lakshminarayana  (1990),  Kiitley  et  al.  (1990),  Leylek  and  Wisler 
(1990),  Moore  and  Moore  (1990),  Subramaniam  (1989),  Warfield  and  Lakshminarayana 
(1989),  for  example]. 

The  comparative  performance  of  implicit  and  explicit  numerical  schemes,  in  three- 
dimensional  steady  flow  computations  is  not  definitive.  Both  explicit  and  implicit 
procedures  generally  incorporate  local  timestepping,  where  estimated  optimum  timesteps 
are  used.  For  explicit  schemes,  this  corresponds  to  operational  Courant-Friedrichs-Lewy 
(1928),  or  CFL,  numbers  just  below  that  allowed  for  numerical  stability  (=  2.8  for 
standard  four-stage  Runge-Kutta  schemes).  In  multidimensional  problems,  implicit  time- 
marching  schemes  are  generally  factored  in  nature,  and  factorization  gives  rise  to  optimum 
local  CFL  numbers  of  the  same  order  (in  three-dimensions,  typical  optimum  CFL  numbers 
of  approximately  5  have  been  reported  for  LU  schemes  [Yokota  and  Caughey  (1987)]  and 
approximately  10  for  ADI  schemes  [Meikle  (1988),  for  example]).  Since  local  timesteps 
are  monotonically  related  to  local  distances  between  grid  points,  the  performance  of  both 
implicit  and  explicit  schemes  suffer  on  the  highly  clustered  grids  required  for  high 
Reynolds  number  flow  computations. 

There  are  several  other  considerations  involved  in  the  choice  between  implicit  and 
explicit  schemes,  including  the  availability  of  steady  state  acceleration  schemes.  Multigrid 
procedures  have  been  used  for  explicit  Euler  calculations  with  great  success  (Ni  (1982),  for 
example)  and  with  some  success  in  implicit  computations  [Yokota  and  Caughey  (1987),  for 
example].  Implicit  residual  smoothing  (IRS)  schemes  can  provide  approximately  a  factor 
of  two  in  CPU  times  savings  in  explicit  procedures,  but  the  performance  of  both  multigrid 
and  IRS  schemes  are  very  sensitive  to  the  grid  stretching  required  in  viscous  flow 
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computations[see  Radespiel  ct  al.  (1990)]. 

Explicit  algorithms,  and  some  implicit  algorithms,  are  fully  vectorizable,  giving 
rise  to  rapid  (CPU  seconds)/(grid  point  *  iteration)  performance  on  advanced  vector 
computer  architectures.  Explicit  procedures  can  be  more  easily  rendered  close  to  100  % 
parallelizable,  so  application  to  the  massively  parallel  architectures,  which  are  now 
emerging,  is  appropriate.  Modified  implicit  solvers,  though,  have  recently  been  shown 
[Wake  and  Egolf  (1991))],  to  provide  efficient  execution  rates  on  such  machines,  as  well. 

Lastly,  the  nature  of  the  discrete  equation  set  being  solved  should  be  considered  in 
choosing  a  particular  numerical  scheme  (or  a  temporal  discretization  procedure  for  certain 
terms).  Specifically,  implicit  treatment  may  improve  damping  properties  when  applied  to 
discrete  transport  equations  with  source  terms.  In  summary,  it  is  currently  not  clear 
whether  implicit  or  explicit  procedures  are  uniformly  more  efficient  in  viscous  three- 
dimensional  turbomachinery  computations. 

An  explicit  numerical  procedure  of  the  Runge-Kutta  class  has  been  chosen  in  this 
thesis.  This  scheme  is  preferred  in  this  work  for  several  reasons.  First,  the  numerical 
formulation  is  somewhat  simpler  than  for  implicit  schemes,  which  is  helpful  from  a  code 
development  standpoint,  when  complications  due  to  additional  physical  modelling  and 
application  to  turbomachinery  configurations  is  required.  Second,  Runge-Kutta  schemes 
are  more  efficient  [Jameson  et  al.  (1981)]  and  naturally  less  dissipative  than  other  explicit 
schemes  (central  and  upwind  difference;  see  Figure  3.2  and  accompanying  discussion), 
allowing  the  user  to  control  more  precisely  levels  of  artificial  dissipation  which  may  corrupt 
solution  accuracy.  Furthermore,  a  number  of  researchers  have  successfully  applied  Runge- 
Kutta  procedures  in  viscous  three-dimensional  turbomachinery  computations  [see  Chima 
and  Yokota  (1988),  Kirtley  et  al.  (1990),  Subramaniam  (1989)  for  example].  Lastly,  as 


6 


will  be  shown,  the  procedures  developed  in  this  thesis  have  been  very  successful,  which  in 
retrospect,  ultimately  justifies  their  choice. 

Due  to  the  wide  range  of  length  scales  associated  with  turbulent  motions  in  high 
Reynolds  number  turbomachinery  flows,  models  which  require  closure  approximations  for 
the  time  averaged  governing  equations  will  be  required  for  computation  of  these  flows  for 
many  years. 

Closure  models  can  be  broadly  classified  into  first  order  and  second  order 
closures.  In  first  order  closures,  the  Reynolds  stresses  (second  order  correlations)  are 
modelled  as  functions  of  gradients  of  the  mean  velocity  field.  Models  which  fall  under  this 
classification  are  commonly  referred  to  as  eddy  viscosity  models,  and  include  algebraic  and 
one-equation  mixing  length  models,  and  two-equation  models.  In  second  order  closures, 
the  transport  equations  for  the  Reynolds  stresses  themselves  are  treated,  where  the  third 
order  correlations  are  expressed  in  terms  of  mean  flow  quantities  and  the  Reynolds  stresses 
themselves.  Models  which  fall  under  this  classification  include  algebraic  Reynolds  stress 
models  (ARSM)  and  full  Reynolds  stress  models  (FRSM). 

The  first  and  simplest  of  turbulence  models  is  the  mixing  length  model  of  Prandtl 
(1925).  This  model  assumes  that  characteristic  time  scales  of  turbulent  transport  by  eddy 
motion  are  of  the  same  order  as  time  scales  of  the  local  mean  flow.  In  simple  shear  layers, 
where  only  one  mean  strain  rate  and  Reynolds  shear  stress  component  are  dominant,  the 
model  takes  the  following  form  : 


-pu'iu2  =  p, 


9Ui 

9x2 


=  CpL 


20Ui 

9x2 


[1.1] 
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where  an  eddy  viscosity  has  been  introduced,  following  Boussinesq  (1877),  who  proposed 
that  the  turbulent  shear  stresses  can  be  modelled  as  being  linearly  related  to  mean  strain  rate 
as  molecular  stresses  are  in  a  Newtonian  fluid.  In  equation  [1.1],  L,  represents  a  local 
"mixing  length,"  which  must  be  prescribed  by  a  combination  of  analytical  and  empirical 
results,  and  C  is  an  empirical  constant.  The  accuracy  of  mixing  length  models  depend  on 
the  accuracy  of  the  mixing  length  prescription.  In  simple  two-dimensional  shear  layers, 
local  mixing  length  distributions  are  well  established,  but  the  evaluation  of  L  becomes 
difficult  in  non-simple  shear  layers,  most  especially  in  internal  flows  such  as 
turbomachinery.  Also,  the  mean  strain  only  provides  a  good  representation  of  the 
turbulence  time  scale  when  the  Reynolds  stress  production  and  dissipation  are 
approximately  in  balance.  Additionally,  transport  and  history  effects  on  local  turbulence 
structure  are  neglected  in  mixing  length  models.  Regardless  of  these  shortcomings, 
extensions  of  equation  [1.1]  to  multidimensional  problems  have  been  very  popular.  In 
such  approaches,  a  tensor  form  extension  of  tire  Boussinesq  approximation  is  employed  : 

-pu'iuj  =  2^itSij  -  |pk5ij  [12] 

where  Sjj  is  the  mean  strain  rate  tensor,  and  k  is  the  turbulent  kinetic  energy.  The  very 
popular  two-layer  model  of  Baldwin  and  Lomax  (1978)  falls  in  this  class,  where  the  last 
term  in  equation  [1.2]  is  neglected,  and  the  eddy  viscosity  is  algebraically  related  to  the 
mixing  length  in  a  fashion  analogous  to  equation  [1.1]. 

If  the  turbulence  kinetic  energy,  k,  is  known  locally,  then  time  scale  arguments 
similar  to  those  justifying  the  mixing  length  hypothesis  yield  : 


Ht  =  Cnpk1/2L 


[1.3] 
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where  now  k1/2  is  taken  as  a  velocity  scale  representative  of  the  turbulence  eddy  transport. 
In  one  equation  models,  a  differential  transport  equation  is  solved  for  k,  but  the  mixing 
length  is  still  prescribed  algebraically.  These  models  have  been  popular  in  external  flow 
computations,  where  mixing  length  prescriptions  are  not  as  ambiguous  as  in  internal  flows 
(King  (1987),  for  example). 


It  was  mentioned  that  the  specification  of  a  local  mixing  length  is  difficult  in 
turbomachinery  flows  where  shear  layers  can  be  highly  three  dimensional,  local  flow 
reversal  can  exist  (in  which  case  equation  [1.1]  provides  that  *uiu2  0  at  locations  within 

the  separation  zone  where  "x2  ),  and  solid  boundary  proximity  can  be  ambiguous. 

For  these  reasons,  two-equation  models,  which  solve  conceptually  a  transport  equation  for 

the  turbulence  length  scale,  have  been  used  extensively.  A  variety  of  two-equation  models 

exist,  all  of  them  have  in  common  that  the  characteristic  local  turbulence  velocity  and  length 

scales,  k1^2,  L,  can  be  extracted  from  the  two  transport  variables.  The  most  well  known  of 

these  is  the  k-e  model,  which  is  based  on  the  pioneering  work  of  a  number  of  researchers 

including  Jones  (1971)  and  Launder  and  Spalding  (1974).  The  underlying  assumption  of 

R ,  -VkL 

i _ i  _ i _ i _ n _ u _ i _ I  .. 


these  models,  is  that,  at  high  enough  local  turbulence  Reynolds  number. 


v  ,  the 


rate  of  change  of  local  turbulence  kinetic  energy  is  equal  to  the  energy  dissipated  in  the 
small  scales  (local  equilibrium  assumption).  Since  the  dynamics  of  the  large  scale 
turbulence  motions  can  be  characterized  by  velocity  and  length  scales,  k1^  and  L, 


dimensional  analysis  provides  that 


e  =  E(k,L)  =  C  ^ 

L 


3U;  3U; 


where  the  dissipation  rate  is  defined  here  as  ^xj  ^xj,  which  is  the  homogeneous 
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component  of  the  total  dissipation  rate  (see  discussion  in  Chapter  2).  Accordingly,  we 
L  =  C^ 

have  e  ,  and  invoking  equation  [1.3],  the  Prandd-Kolmogorov  relation  for  eddy 

viscosity  becomes : 


[1-5] 


where  is  an  empirical  constant. 

Launder  and  Spalding  (1974)  provided  the  so-called  standard  high  Reynolds 
number  k-e  model  and  recommended  empirically  obtained  values  for  the  modelling 
constants  which  appear.  In  order  to  account  for  near  wall  effects,  where  the  high  Reynolds 
number  form  of  the  k  and  e  equations  are  not  valid,  Launder  and  Spalding  (1974) 
introduced  wall  functions  which  incorporate  assumed  near  wall  log-law  velocity  profiles. 
Accordingly,  wall  functions  are  applied  at  the  first  grid  point  off  of  the  wall  which,  for 
consistency,  should  lie  in  the  inertial  sublayer  (nominally  30  <  y+2  <  300).  Boundary 
conditions  on  k  and  e  which  assume  local  equilibrium  (production  of  k  =  dissipation)  are 
also  applied  at  this  location.  Wall  functions  have  the  advantage  that  significant  computer 
resource  savings  can  be  realized  since  near  wall  physics  need  not  be  resolved  explicitly. 
This  has  kept  wall  functions  popular,  even  in  state-of-the-art  three-dimensional 
turbomachinery  flow  solvers  (Leylek  and  Wisler  (1990),  Bansod  and  Rhie  (1990),  for 
example). 

However,  since  wall  functions  rely  on  the  assumption  of  a  locally  logarithmic 
velocity  distribution  and  the  local  equilibrium  of  the  turbulence,  their  accurate  applicability 
is  severely  restricted.  Specifically,  in  locally  separated  flows,  no  log-law  region  exists, 
and  in  three-dimensional  flows,  secondary  motions  (say  spanwise  flows)  often  extend  well 
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into  the  inner  layer,  giving  rise  to  skewed  boundary  layer  profiles.  Additionally,  implicit  in 
the  use  of  wall  functions  is  that  blade  boundary  layers  are  fully  turbulent  Lastly, 
application  of  wall  functions  in  turbomachinery  configurations  gives  rise  to  several 
ambiguities,  including  leading  edge  and  blade-endwall  intersection  treatment  and  wake 
centerline  treatment.  For  these  reasons,  the  wall  function  approach  is  not  adopted  in  the 
present  work. 

Low  Reynolds  number  forms  of  the  k-e  model  directly  model  near-wall  influences 
on  the  turbulence  structure  by  modifying  the  k  and  e  equations.  Specifically,  they  are 
distinguished  from  high  Reynolds  number  forms  in  that  several  modelling  parameters 
which  account  for  the  proximity  of  solid  boundaries  are  accounted  for  (see  Patel  et  al. 
(1985)  for  example).  Specifically,  molecular  viscosity  is  retained  in  the  diffusion  terms, 
model  constants  in  the  dissipation  rate  equation  are  factored  with  near  wall  modification 
functions,  and  in  some  cases  near  wall  source  functions  are  incorporated  to  model  directly 
the  behavior  of  the  dissipation  rate  near  the  wall  (refer  to  equations  [2.10]  -  [2.12]).  The 
eddy  viscosity,  and  thereby  conceptually  the  Reynolds  stresses,  are  also  damped  as  the 
wall  is  approached  in  a  fashion  consistent  with  Van  Driest  (1956)  damping.  Specifically, 
the  Prandtl-Kolmogorov  relation  is  modified  as  : 

=  f|iC|iP^  [1.6] 

where  f^  provides  appropriate  exponential  damping  as  the  wall  is  approached. 

Low  Reynolds  number  models  have  the  advantage  that  near  wall  boundary  layer 
skew  due  to  secondary  motions  can  be  resolved.  This  is  important  in  flows  such  as 
centrifugal  compressors,  where  significant  cross  flow  velocity  components  can  exist  very 


11 


close  to  the  wall  (refer  to  results  in  Chapter  6).  These  models  can  be  applied  in  an 
unambiguous  fashion  in  complex  geometries  (all  that  is  needed  is  an  appropriate  wall 
proximity  length  scale  prescription;  see  Chapter  2).  This  is  important  especially  at 
turbomachinery  boundaries  such  as  the  interface  between  rotating  and  stationary  portions  of 
a  hub,  a  stationary  blade  adjacent  to  a  rotating  (in  the  relative  frame)  shroud,  and  the 
interface  between  solid  and  periodic  boundaries,  such  as  at  leading  and  trailing  edges  and 
tip  clearance  gaps.  In  the  opinion  of  the  author,  the  use  of  wall  functions  along  these 
boundaries  would  be  difficult  and  inappropriate.  Lastly,  low  Reynolds  number  k-e  models 
can  model  the  influence  of  ffeestream  turbulence  on  turbulence  transition  (Jones  and 
Launder  (1972)  for  example),  albeit  only  if  the  modelled  flow  is  such  that  transition  is 
strongly  influenced  by  local  ffeestream  intensity. 

The  improvements  in  prediction  accuracy  which  can  be  obtained  using  two- 
equation  models  over  mixing  length  models  is  well  documented.  For  example,  the  k-e 
model  has  been  shown  to  provide  better  predictions  than  algebraic  models  for  flows  with 
adverse  pressure  gradient  [Kirtley  and  Lakshminarayana  (1985)]  and  for  2-D  shock- 
boundary  layer  flow  on  curved  surfaces  [Degrez  and  VanDromme  (1985)].  However, 
there  are  still  well  acknowledged  deficiencies  [see  Lakshminarayana  (1986]  and  Speziale 
(1989),  for  example)  associated  with  the  use  of  two-equation  models.  Such  models  are 
eddy  viscosity  models,  and  as  such  invoke  the  Boussinesq  approximation  which  implies 
that  the  Reynolds  stress  tensor  is  aligned  with  the  mean  strain  rate  tensor.  Accordingly, 
eddy  viscosity  models  cannot  predict  secondary  flows  of  the  second  kind,  which  are 
induced  by  Reynolds  intensity  anisotropies,  and  are  present  even  in  unidirectional  flows 
(see  further  discussion  in  Chapter  5,  Section  1).  Also,  eddy  viscosity  models  cannot 
account  for  component  Reynolds  stress  relaxation/amplification  effects;  instead  they  predict 
instantaneous  return  to  isotropy  in  the  absence  of  local  mean  strain.  Another  criticism  of 
these  models,  is  that  unless  empirically  modified,  they  cannot  account  for  the  influence  of 
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local  streamline  curvature  and  system  rotation  on  local  turbulence  structure. 

Full  Reynolds  stress  models  (FRSM)  of  the  form  proposed  by  Launder  et  al. 
(1975) : 


Du 


L  a  .  I  —du,  —  dui 

D^  =  iCsa^EUkU>-a^)| +  (-  UiUk5S^  -  “JUk5^ 


+  {'^-l^(UiUj  '  3^>Jk)  *^ij  *  +  ^ij.w)  +  {  ‘  3^>je 


[1.7] 


or : 


Convection  =  {Diffusion,  Djj]  +  {Production  by  mean  strain  rate,  P,j) 

+  {Pressure  strain)  +  {Dissipation} 

can  in  principle  account  for  the  eddy  viscosity  modelling  deficiencies  alluded  to  above.  In 
these  closures,  the  modelled  transport  equations  for  the  six  individual  Reynolds  stress 
components,  equation  [1.7],  are  solved,  along  with  an  equation  for  turbulence  energy 
dissipation  rate.  Various  FRSM  have  been  shown  to  account  correctly  for  Reynolds  stress 
anisotropy  phenomena  and  amplification/relaxation  effects,  and  account  intrinsically  for 
curvature  and  rotation  effects  since  fluctuating  momentum  balance  information  is 
incorporated  in  the  exact  equations  (see  Gibson  and  Rodi  (1978)  and  Launder  et  al.  (1987) 
for  example).  Currently,  however,  the  use  of  FRS  closures  has  been  limited  to  simple 
flows  wherein  reduced  forms  of  the  models  and  less  computationally  intensive  numerical 
schemes  are  applicable.  This  is  due  to  the  computational  intensity  of  solving  seven 
transport  equations,  each  of  which  contain  complicated  source  terms,  in  addition  to  the 
mean  flow  equations. 


Rodi  (1976)  proposed  an  approximation  which  reduces  the  modelled 
incompressible  full  Reynolds  stress  (FRS)  equations  to  algebraic  form.  Specifically,  he 
presented  the  FRSM  of  Launder  et  al.  (1975),  equation  [1.7],  and  proposed  the  following 
approximation : 


where  Dk  represents  the  viscous  diffusion  term  in  the  turbulence  kinetic  energy  equation 
(1/2  the  trace  of  equation  [1.7]).  This  approximation  assumes  that  the  convective  and 
diffusive  transport  of  individual  Reynolds  stress  components  are  locally  proportional  to 
transport  of  turbulent  kinetic  energy,  the  constant  of  proportionality  being  the  local  ratio  of 
the  Reynolds  stress  component  to  the  turbulent  kinetic  energy.  This  simplification  reduces 
the  system  of  six  non-linear  Reynolds  stress  partial  differential  equations  to  a  system  of  six 
non-linear  algebraic  equations  in  the  six  unknown  stresses.  The  advantage  to  such  an 
approach  is  that  solution  of  the  algebraic  systems  should  be  less  computationally  intensive, 
though  solution  of  a  6  x  6  system  of  non-linear  algebraic  equations  at  every  grid  point,  at 
every  iteration,  can  be  computationally  intensive.  Also,  transport  equations  for  k  and  e 
must  be  solved  in  order  to  provide  the  values  which  appear  in  equation  [1.8]. 

In  this  thesis,  a  low  Reynolds  number  form  of  the  k-e  model  is  used.  The 
particular  model  chosen  is  a  three-dimensional,  compressible  extension  to  a  model 
proposed  by  Chien  (1982).  This  model  was  one  of  four  low  Reynolds  number  models 
deemed  to  perform  well  in  the  prediction  of  a  variety  shear  layer  flows,  in  a  recent  review 
article  by  Patel  et  al.  (1985).  The  formulation  for  this  model  is  presented  in  equations 
[2.10]  -  [2.12].  Also,  a  hybrid  low  Reynolds  number  k-e/high  Reynolds  number  ARS 
model  has  been  developed  and  applied.  One  of  the  goals  of  this  thesis  work  has  been  to 
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compute  the  flow  in  a  centrifugal  compressor.  It  was  felt  that  for  adequate  resolution  of 
the  strong  secondary  motions  in  such  machines,  a  low  Reynolds  number  k-e  model  would 
be  significantly  more  appropriate  than  a  mixing  length  model  or  a  high  Reynolds  number 
form  with  wall  functions.  Application  of  the  model  adopted  to  the  complex  geometries  in 
impellers,  and  all  turbomachinery  applications,  is  relatively  straightforward.  Additionally, 
as  second  order  closures  are  adopted  in  the  code,  transport  equations  for  k  and/or  e  need  to 
be  solved  anyway,  so  the  implementation,  stability  and  application  results  presented  in  this 
thesis  are  directly  relevant  to  incorporation  of  ARS  and  FRS  closures. 

The  flow  in  centrifugal  compressors  and  pumps  provide  one  of  the  most 
challenging  problems  available  to  internal  flow  aerodynamicists.  The  complex  three- 
dimensional  physics  in  all  of  these  machines  has  defied  accurate  prediction  to  date.  The 
real  flow  in  an  impeller  is  drastically  different  from  the  desirable  inviscid  flow  arising  from 
the  influence  of  curvature  and  rotation  on  rotor  streamtubes  in  the  absence  of  shear.  A 
dominant  boundary  layer  separation,  or  "wake"  appears  in  many  radial  and  mixed  flow 
centrifugals,  which  alters  strikingly  the  flow  structure  in  the  radial  section  of  the  impeller. 
A  variety  of  strong  secondary  flows  develop  and  interact  due  to  the  influence  of  centrifugal 
and  coriolis  forces  on  blade  and  endwall  shear  layers.  Curvature  and  rotation  also 
influence  local  turbulence  intensity  levels,  thereby  affecting  the  structure  of  developing 
boundary  layers.  Tip  clearance  flows  and  unsteady  impeller-diffuser  interactions  introduce 
additional  anomalies  to  an  already  complex  passage  flow.  Additionally,  there  can  be 
compressibility  effects,  even  shock  surfaces  in  the  passage.  Regions  of  flow  reversal  can 
be  present  at  off  design  operation.  In  short,  the  complete  physics  involved  are  highly 
complex,  and  therefore,  Navier-Stokes  analyses  to  date  have  provided  only  qualitative 
accuracy. 


Full  Navier-Stokes  procedures  have  been  implemented  to  compute  centrifugal 


compressor  flows  by  several  researchers  recently  (refer  to  Dawes  (1988),  Choi  and  Knight 
(1989),  Goto  (1990),  Bansod  and  Rhie  (1990),  Hah  and  Krain  (1989),  Domey  and  Davis 
(1990)).  Most  of  these  researchers  have  extracted  qualitative  comparison  with 
experimentally  obtained  pressure  distributions,  and  typically  provide  reasonable 
performance  parameter  trends  (with  varying  mass  flow  rate  or  tip  clearance  gap  for 
instance).  Very  popular  is  the  use  of  solution  contour  plots  of  meridional  velocity  and 
relative  stagnation  pressure,  and  cross  flow  velocity  vector  plots.  Such  numerical 
simulation  results  provide  insight  into  the  variety  of  complex  flow  phenomena  which  exist 
in  these  machines,  and  substantiate  the  ability  of  Navier-Stokes  procedures  to  qualitatively 
model  these  physics.  However,  direct  comparison  of  predicted  mean  internal  passage  flow 
properties  to  measured  results  are  consistently  not  provided,  since  quantitatively, 
predictions  are  generally  poor.  Additionally,  only  two  of  these  works  (Choi  and  Knight 
(1989),  Hah  and  Krain  (1989))  have  adopted  the  degree  of  turbulence  modelling 
complexity  of  a  low  Reynolds  number  two-equation  model. 

In  the  last  chapter  of  this  thesis,  the  results  and  interpretation  of  a  transonic 
impeller  computation  is  provided.  Direct  comparisons  of  predicted  meridional  velocities  are 
made  with  L2F  measurements  due  to  Krain  (1988).  Relative  helicity  density  is  introduced 
as  an  excellent  parameter  for  interpretation  of  secondary  flow  development  Extensive 
physical  interpretation  is  provided  for  a  solution  in  which  the  low  Reynolds  number  k-e 
model  is  used.  Additionally,  results  using  a  hybrid  low  Reynolds  number  k-e/ARS  model 
are  presented  which  indicate  that  although  meridional  mean  flow  is  not  significantly 
influenced  by  the  modelling  choice,  cross  flows,  wall  shear  stress  and  machine 
performance  parameters  are. 

The  choice  of  an  explicit  numerical  procedure  in  conjunction  with  a  low  Reynolds 
number  k-e  model  is  not  without  its  consequences.  Specifically,  nearly  as  well 
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documented  as  the  accuracy  advantages  of  using  low  Reynolds  number  form  k-e  models 
over  algebraic  eddy  viscosity  models,  are  the  perceived  numerical  disadvantages,  and 
thereby  the  unpopularity,  of  implementing  two-equation  models  in  time-marching  Navier- 
Stokes  procedures  (Coakley  (1983),  for  example).  For  these  reasons,  algebraic  eddy 
viscosity  models  are  often  used  to  approximate  the  apparent  stresses  in  CFD  codes,  since 
these  models  have  little  computational  overhead  and  do  not  adversely  affect  the  stability  of 
the  scheme.  When  transport  models  are  used,  the  numerical  stiffness  associated  with 
implicit  and  explicit  computation  of  the  Navier-Stokes  and  turbulence  transport  equations 
on  the  highly  stretched  grids  required  by  low  Reynolds  number  two-equation  models  has 
prompted  researchers  to  adopt  a  variety  of  strategies  to  obtain  convergent  solutions.  Wall 
functions  (Eliasson  (1988),  Grasso  and  Speziale  (1989),  Holmes  and  Connell  (1989)  for 
instance),  require  less  near  wall  grid  clustering  and  bridge  much  of  the  inner  region  of  the 
boundary  layer,  where  turbulence  source  terms  are  large.  Liu  (1987)  has  used  a  hybrid 
approach,  where  an  algebraic  Van-Driest  damping  model  is  used  to  model  the  near-wall 
physics,  thereby  also  relieving  these  two  stability  constraints.  Pointwise  implicit  source 
term  treatment  has  been  used  in  implicit  (Coakley  (1983),  Sahu  and  Danberg  (1986), 
Martelli  and  Michelassi  (1990))  and  explicit  (Liu  (1987),  Gerolymos  (1990),  Mavriplis  and 
Martinelli  (1991))  codes.  Common  to  all  of  these  approaches  is  their  success  in  predicting 
a  variety  of  complex  turbulent  flows  in  convergent  fashion.  In  addition,  most  of  these 
researchers  cite  the  need  for  ad  hoc  stabilization  procedures  in  the  early  stages  of  iteration 
where  source  terms  can  be  very  large  if  the  flow  is  initialized  in  a  simple  fashion. 

The  above  considerations  have  motivated  the  comprehensive  stability  investigation 
provided  in  Chapter  3  of  this  thesis,  in  which  the  stability  of  the  turbulence  model  is  one  of 
the  emphasized  items.  The  author  is  familiar  with  two  attempts  in  the  literature  to  analyze 
the  stability  of  the  discrete  k-e  equations.  Liu  (1987),  analyzed  the  2D  "uncoupled"  (see 
terminology  discussion  in  Chapter  3),  two-equation  k-e  system  and  provided  an  expression 
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for  a  local  timestep  based  on  implicit  source  term  treatment  using  a  2-stage  Runge-Kutta 
procedure.  Eliasson  (1988)  provided  a  VonNeumann  analysis  of  the  coupled  five  equation 
incompressible  2D  Navier-Stokes/k-e  system.  In  his  work,  the  influence  of  turbulence 
production  on  the  stability  of  the  scheme  was  neglected.  Eliasson  thereby  obtained  the 
result  that  the  remaining  dissipation  terms  effectively  provide  damping  to  the  scheme.  He 
reconciles  this  result  with  the  observation  that  it  is  not  necessary  to  add  fourth  order 
artificial  dissipation  to  the  k  or  e  equations,  which  has  also  been  the  experience  of  this 
author.  Though  this  neglect  of  production  with  respect  to  dissipation  in  the  stability 
analysis  may  be  useful  when  using  wall  functions,  as  Eliasson  did,  it  was  felt  that  it  would 
not  be  appropriate  to  adopt  this  simplification  in  the  present  work.  In  Chapter  3,  much 
more  general  stability  results  are  obtained  than  by  these  authors,  and  numerical  verification 
studies  are  performed  to  corroborate  these  analyses. 

The  explicit  flow  solver  used  for  the  numerical  studies  included  in  this  thesis  was 
written  and  developed  by  the  author  and  is  described  briefly  here.  The  code  is  a  structured 
H-grid,  explicit,  compressible,  3D,  full  Navier-Stokes  flow  solver  which  incorporates  a 
compressible  low  Reynolds  number  k-e  model.  Inlet  characteristic  boundary  conditions, 
system  rotation  and  periodic  boundaries  are  available  to  accommodate  turbomachinery 
applications.  Euler  and  periodic  boundaries  consistent  with  either  rectilinear  or  annular 
geometries  are  available.  The  modified  eigenvalue  and  velocity  artificial  dissipation 
scalings  presented  in  Chapter  3  are  implemented.  An  algebraic  Reynolds  stress  model 
which  may  be  hybrid  with  the  low  Reynolds  number  k-e  model  is  available  in  the  code. 
Code  performance  issues  are  discussed  in  Chapter  3. 

In  Chapter  2,  the  governing  equations  for  the  mean  flow  and  turbulence  models  are 
provided.  Chapter  3  provides  details  of  the  spatial  and  temporal  discretization  procedures 
used.  This  includes  the  introduction  of  a  compact  central  finite  difference  flux  evaluation 
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scheme  which  eliminates  the  metric  singularity  at  the  interface  between  solid  and  periodic 
boundaries,  often  attributed  to  finite  difference  procedures.  Also  included  in  Chapter  3  are 
extensive  analyses  and  discussion  of  artificial  dissipation  issues  which  arise  in  the  high 
Reynolds  number  viscous  flow  computations  of  interest  in  this  thesis.  The  motivation  and 
stability  implications  of  numerically  coupling  the  turbulence  transport  equations  to  the  mean 
flow  equations  are  then  analyzed  and  discussed.  Next,  a  vector  VonNeumann  analysis  of 
the  seven  discrete,  numerically  coupled  governing  equations  is  undertaken,  and  some  local 
stability  results  obtained.  Order  of  magnitude  arguments  and  numerical  verification 
experiments  are  applied  to  these  results  from  which  some  conclusions  of  relevance  to  the 
engineering  analyst  are  drawn.  ARSM  numerical  implementation  and  stability  results,  as 
well  as  overall  code  performance  issues  are  also  presented  in  Chapter  3. 

Chapter  4  provides  turbulence  modelling  validation  and  two-dimensional  cascade 
computation0.  Results  and  interpretation  are  included  for  a  supersonic  compressor  cascade 
operating  at  unique  incidence,  and  a  low  subsonic  double-circular-arc  cascade  with  mean 
flow  reversal.  The  flowfields  within  each  of  these  have  been  investigated  experimentally, 
and  these  results  are  used  for  comparison.  In  Chapter  5,  the  three-dimensional 
methodology  is  applied  to  a  turbulent  duct  flow  for  which  CFD  validation-quality 
experimental  data  is  available.  Also,  results  of  a  large  scale  three-dimensional  compressor 
rotor  flow  computation  are  presented.  This  Navier-Stokes  solution  is  compared  with  a 
variety  of  rotor  and  wake  flowfield  measurements,  and  extensive  interpretation  is  provided. 
As  mentioned  above,  Chapter  6  details  the  computation  of  a  centrifugal  compressor 
flowfield  using  the  present  method. 
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Chapter  7  summarizes  the  conclusions  arrived  at  in  this  thesis.  The  range  and 
applicability  of  the  methodologies  presented  are  also  discussed  there,  as  are  limitations  and 
suggestions  for  further  study. 
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CHAPTER  2 

THEORETICAL  FORMULATION 

2.1  FavTe  Averaged  Governing  Equations 


In  the  present  development,  the  density-weighted  time-averaging  developed  by 
Favre  (1965)  is  used.  This  decomposition  has  advantages  in  flow  computations  with 
variable  density  (see  Jones  (1980)).  Specifically,  the  averaged  governing  equations  are  of 
simpler  form  and  the  physical  interpretation  of  terms  in  the  equations  is  clearer  than  when 
conventional  time  averaging  is  used.  Reynolds  (time)  averaging,  defined  for  a  scalar,  <J>,  as 


<j>  =  <{>  +  <J>'  ,  where  0  = 


lim  i 

t— »oo 


tf 


0  ds 


[2.1] 


is  used  for  pressure,  density,  molecular  stress  tensor  and  molecular  heat  flux  vector.  Favre 
(density-weighted  time)  averaging,  defined  for  scalar,  0,  as 


0  =  0  +  0"  ,  where  0  =  ^ 

P  ,  [2.2] 

is  used  for  velocity  components,  internal  energy,  turbulent  kinetic  energy  and  turbulent 
energy  dissipation  rate. 


2.1.1  Mean  Flow  Equations 

Applying  the  density-weighted  time-averaging  procedure  to  the  continuity, 
momentum  and  energy  equations  allows  one  to  write  the  five  mean  flow  equations  in 
conservation  form,  in  cartesian  tensor  notation,  as : 
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where  the  effective  stress  tensor  and  effective  heat  flux  vector  are  given  by : 


[2.3] 
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[2.4] 


[2.5] 


All  velocity  components  in  equations  [2.3]  -  [2.5]  are  relative  to  a  reference  frame 
rotating  with  constant  angular  velocity,  coj  (s1).  In  the  following  developments,  co  is 
taken,  without  loss  of  generality,  to  be  coincident  with  the  x-axis.  The  particular  choice  of 
energy  variable,  P®0R  =  P®  +  pW  12  -  poA2^  ^as  jj,e  advantage  that  for  the  inviscid  only 
problem,  the  steady  state  converged  energy  equation  reduces  to  a  statement  of  conservation 
of  rothalpy,  I,  along  streamlines.  The  apparent  body  force  terms,  which  arise  in  the 
momentum  equations  due  to  system  rotation,  are  linear  in  relative  velocity,  and  therefore  do 
not  introduce  additional  stresses  due  to  the  density- weighted  time-averaging  procedure.  In 
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the  derivation  of  equation  [2.3],  it  has  been  assumed  that  the  time  averaged  molecular  stress 
tensor  and  the  molecular  heat  flux  vector  are  equivalent  to  their  density-weighted  time 
averaged  values.  This  approximation  should  be  a  good  one  since,  compared  to  turbulent 
diffusion,  molecular  diffusion  is  only  significant  near  solid  boundaries,  where  local  Mach 


number  and  hence  density  fluctuations  are  small.  For  the  compressible  flows  considered 
in  this  thesis,  a  perfect  gas  equation  of  state  is  provided  as  P  =  where 

T  =  -LkR-2i  +  S!£!i) 

Cv\  2  2  /.  It  is  further  noted  that  by  definition  of  the  energy  transport 


variable  used. 


2  ,  Favre  averaging  gives  rise  to  a  turbulence  energy  term. 


-  -  Cyf  ,  W2  ,  '  coV 

ie.  v  2  “  2  ,  where  turbulent  kinetic  energy  is  defined  as 


~  ipui'ii- 

k  =  — -  -rr - 


Strictly,  this  term  should  be  accomodated  in  pressure  and  temperature  terms  in  the  inviscid 
flux  vectors  in  the  momentum  equation  and  in  the  viscous  flux  vectors  in  the  energy 
equation.  These  terms  are  usually  neglected  in  the  literature  and  in  this  thesis,  though 
Wilcox  (1991)  has  pointed  out  that  these  terms  become  important  above  local  Mach 
numbers  near  3. 


2.1.2  Eddv  Diffusivitv  Approximations 


Equations  [2.3]  contain  unknown  correlated  fluctuation  quantities  which  arise 
from  averaging  the  nonlinear  convection  operators  in  the  momentum  and  energy  equations. 
Eddy  viscosity  models  employ  a  Boussinesq  (1877)  approximation,  where  the  Reynolds 

stress  tensor,  Puiuj  j  is  assumed  to  be  linearly  proportional  to  the  strain  rate  tensor  : 


■pu'u'  =  + 

1  J  |]9xj  dxj)  3  J9xk 


fVk 


[2.6] 


where  an  eddy  viscosity,  pt,  has  been  introduced.  This  is  one  of  the  approaches  taken  in 
the  present  work. 
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Analogous  to  assuming  a  linear  relationship  between  apparent  stress  and  velocity 
gradient,  a  gradient  diffusion  hypothesis  is  adopted  which  assumes  a  linear  relationship 

"  0 

between  apparent  heat  flux,  Puie  ,  and  internal  energy  gradient : 


'  *  de 

Puje  -  ‘  Y  Pt  energy^ 

or  introducing  a  turbulent  Prandtl  number,  Prt  =  (it/(j.t  energy* 


0  0 


pu;e  =  -  y 


Pt8e 

Prt9xi 


[2.7] 


The  effective  stress  tensor  and  effective  heat  flux  vector  provided  in  equation  [2.4]  can  now 
be  written : 


where  e  has  been  replaced  by  CvT  in  the  present  work,  air  is  used  as  the  working  fluid  in 
all  test  cases.  Accordingly,  the  laminar  Prandtl  number,  Prj,  is  set  to  0.72.  The  turbulent 
Prandtl  number,  Prt,  is  set  to  a  value  of  0.91  (Cebici  (1970)  for  example). 


To  close  the  system  of  governing  equations,  [2.3],  it  remains  to  obtain  the  spatial 
distribution  of  eddy  viscosity  and  turbulent  kinetic  energy. 
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In  low  Reynolds  number  k-e  turbulence  models,  the  eddy  viscosity  is  obtained 


from  the  Prandtl-Kolmogorov  relation : 


CjifjiP  k 


~  tPUiUj 
k  s  2.  _  . 


where  the  turbulent  kinetic  energy  is  defined  as  P  ,  the  homogeneous  component 


9u;  9u- 

vp3^ 


of  the  turbulence  energy  dissipation  rate  is  defined  as  P  ,  =  0.09  (obtained 

empirically  (Launder  ans  Spalding  (1974)))  and  f^  is  a  damping  function  to  accommodate 
the  direct  influence  of  solid  boundary  proximity  on  the  local  turbulence  structure. 


In  general,  low-Reynolds  number  forms  of  the  k-e  model  can  b"  Favre  averaged 
and  cast  in  conservation  form  for  compressible  flows  as  (see  Coakley  (1983),  Narayan 
(1991)  for  example) : 


¥i+^>=#+^£]+p-pe+i) 

¥+ #+^i]+(f'c'p-f2C^+£ 


[2.10] 


where  production  of  turbulent  kinetic  energy  is  defined  as  : 
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[2.11] 
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and  Prk,  Pre,  flt  f2,  Clt  C2,  D,  and  £  are  modelling  parameters. 

The  particular  form  of  low  Reynolds  number  model  used  in  the  code  was 
originally  devised  by  Chien  (1982)  for  incompressible  flow.  For  this  model,  the  constants 
and  functions  in  equations  [2.9]  and  [2.10]  are  given  by  : 


Prk  =  1.0  ,  Pr£  =  1.3 


f^  =  1  -  exp(-.0U5y+)  ,  fj  =  1  ,  f2  =  1  -  2-exp(-RT/36) 


q,  =  0.09  ,  Cl  =  1.35  ,  C2  =  1.80 


i1  1 1 


RrJX.rJ-f,u--JW 

[tie  W  V  p 


[2.12] 


In  turbomachinery  blade  row  computations,  £  corresponds  to  the  streamwise 
coordinate,  r\  to  the  pitchwise  coordinate,  and  £  to  the  spanwise  coordinate.  The  wall 
proximity  length  scale,  /,  used  in  equations  [2.12],  accounts  for  the  influence  of  adjacent 
walls  (such  as  blade-endwall  comers),  and  is  defined  as 


^  ss  ps  ^  hub  ^  shroud  (  where  and  are  absolute  distances  from  solid 
boundaries  in  the  two  crossflow  directions.  Upstream  and  downstream  of  blade  rows, 

is  replaced  by  absolute  distance  from  leading  or  trailing  edge,  respectively.  In  tip  clearance 
gap  regions,  /^hub  is  replaced  by  / ^  Up. 


In  high  Reynolds  number  k-e  models,  the  dissipation  variable,  e,  is  the  dissipation 

9u-  9u; 

-  vpa^ 

rate  in  homogeneous  turbulence,  P  .  The  transport  variable,  e,  used  in 

Chien's  and  certain  other  low  Reynolds  number  k-e  models  is  the  isotropic  component  of 
the  dissipation.  As  discussed  by  Jones  and  Launder  (1972),  at  high  local  Reynolds 
numbers,  the  non-isotropic  component  of  dissipation  is  negligible,  so  the  model  remains 
valid  in  these  regions.  Near  a  solid  wall,  however,  the  dissipation  is  not  negligible,  though 
the  isotropic  component  of  the  dissipation  goes  to  zero.  The  term,  D,  accounts  for  the  non¬ 
zero  value  of  total  dissipation  near  the  wall,  so  that  the  model  also  remains  valid  near  solid 
boundaries,  but  the  convenience  of  specifying  an  £  =  0  Dirichlet  condition  is  retained. 
Specifically,  then,  D  represents  viscous  diffusion  of  turbulence  kinetic  energy  at  the  wall, 
which  is  in  balance  with  dissipation  there.  Some  numerical  stability  implications  of  this 
choice  of  dissipation  variable  are  presented  in  Chapter  3. 


Since  the  k-E  equations  have  been  cast  in  compressible  form,  the  direct  influence  of 
the  dilitation  terms  in  equations  [2.8]  and  [2.1 1]  are  incorporated.  However,  the  other 
modelling  assumptions  invoked  in  this  thesis  are  essentially  those  for  incompressible  flow. 
Specifically,  terms  in  the  unmodelled  e  equation  which  contain  density  fluctuation  terms, 
p',  are  neglected.  The  unmodelled,  Favre  averaged  e  equation  is  extremely  complicated 
due  to  the  presence  of  many  terms  which  contain  density  fluctuations.  Jones  (1980)  points 
out  that  the  physical  interpretation  of  exact  terms  in  the  incompressible  £  equation  are  not 
helpful  in  modelling  the  equation.  He  argues  that  there  is,  then,  little  point  in  trying  to 
model  the  much  more  complicated  exact  equation,  with  varying  density.  Accordingly, 
most  efforts  to  date  to  model  the  £  equation,  in  compressible  flows,  have  retained  standard 
incompressible  modelling. 
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For  the  density-weighted  time-averaged  turbulent  kinetic  energy  equation,  several 
researchers  have  chosen  to  explicitly  model  some  of  the  extra  terms  which  arise  due  to  the 
Favre  averaging  procedure.  Specifically,  VanDromme  (1983)  and  Grasso  and  Speziale 
(1989)  have  modelled  some  or  all  of  the  extra  pressure  work  correlation  terms,  dilitation, 
and  density  fluctuation  terms  which  arise  in  the  exact  Favre  averaged  k  equation,  to  account 
for  extra  strain  effects  in  mixing  layer  flows  where  buoyancy  effects  are  important,  and  in 
supersonic  adverse  pressure  gradient  flows.  Though  the  framework  for  incorporation  of 
these  types  of  corrections  is  available  in  a  code  using  the  Favre  averaging  procedure,  no 
such  modelling  is  employed  in  this  work. 

Stability  implications  and  validation  of  the  low  Reynolds  number  k-E  model  are 
provided  in  Chapters  3  and  4  respectively. 

Anticipation  of  the  computation  of  a  centrifugal  compressor  flow  (Chapter  6) 
required  that  no  thin  layer  approximations  be  made  in  either  the  mean  flow  or  turbulence 
transport  equations. 

2.1.4  Algebraic  Reynolds  Stress  Model  ('ARSMl 

A  class  of  ARS  models,  [Rodi  (1976),  Galmes  and  Lakshminarayana  (1984), 
Andersson  and  Nilsen  (1989),  for  example]  can  be  derived  by  application  of  equation  [1.8] 
to  modelled  FRS  models  such  as  that  given  in  equation  [1.7],  yielding  the  following 
general  algebraic  expressions  for  the  six  individual  Reynolds  stress  components  : 

-  UjUj  =  -  ^-8ijk  ”  ^Tjj 


where 
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T  _  Rij (Q  -  C3)  +  (Pjj  -  2PSij/3)(l-C2)  +  Oij.w  __  __ 

y  P  +  e(Ci-l)  f  Rik  =  -  2a)p(  eipju’kUj  -  eicpjpu'iuj ) 


=  2P  =  P„ 


[2.13] 


and  the  parameters  Q,  Cj,  C2,  C3,  foj.w  arise  from  the  pressure  strain  correlation 
modelling  in  the  particular  FRSM  chosen.  Ry,  the  production  by  rotation  tensor,  arises 
when  the  Reynolds  stress  transport  equations  are  derived  in  a  rotating  coordinate  system  of 
angular  velocity  vector,  CDp. 


A  number  of  researchers  [Launder  et  al.  (1975),  Galmes  and  Lakshminarayana 
(1984),  Shih  and  Lumley  (1986),  Launder  et  al.  (1987)  for  example]  have  proposed  near 
wall  modelling  for  FRS  or  ARS  closures,  which  account  for  the  influence  of  the  wall  on 
the  pressure  strain  correlation,  ^‘j.w  (the  so  called  "echo  effect,"  which  accounts  for  the 
damping  of  the  wall  normal  velocity  fluctuation  and  redistribution  of  its  energy  in  the  other 
two  intensities),  and/or  very  near  wall  viscous  damping  effects.  In  practical  numerical 
computations,  wall  functions  are  often  used  in  both  FRS  [Shih  and  Lumley  (1986), 
Launder  et  al.  (1987)]  and  ARS  [Warfield  and  Lakshminarayana  (1987),  Zhang  and 
Lakshminarayana  (1990)]  models  to  bridge  the  near  wall  region,  thereby  avoiding  the 
necessity  for  viscous  near  wall  modelling.  Near  wall  pressure- strain  modelling,  though,  is 
retained  in  boundary  layer  code  computations  using  FRSM  by  Shih  and  Lumley  (1986)  and 
Launder  et  al.  (1987),  and  an  ARS  version  of  Shih  and  Lumley's  model  due  to  Zhang  and 
Lakshminarayana  (1990).  For  full-scale  Navier-Stokes  computations  of  the  type 
performed  by  Warfield  and  Lakshminarayana  (1987)  and  of  interest  in  this  thesis,  these 
terms  are  not  retained,  for  two  reasons.  Firstly,  available  models  for  ^>j.w  involve  many 
very  complicated  terms  (Shih  and  Lumley  (1986),  Launder  et  al.  (1987),  for  example).  In 
full  Navier-Stokes  implementation,  the  computation  of  this  term  alone  would  significantly 


burden  the  overall  CPU  performance  of  the  method.  The  second  reason  for  the  neglect  of 
this  near  wall  pressure  strain  term  is  its  implicit  incorporation  in  the  hybrid  modelling 
approach  presented  below. 
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In  this  thesis,  a  hybrid  low  Reynolds  number  k-e  model/high  Reynolds  number 
ARSM  approach  has  been  devised.  In  the  fully  turbulent  region  (y+  =  O(102)),  a 
compressible  extension  to  the  high  Reynolds  number  form  of  the  ARSM  due  to  Galmes 
and  Lakshminarayana  (1984)  is  adopted: 

-  pu-u]  =  -  jSijpk  -  pkTij 

where 


T  _  Rij(2-C2)/2  +  (Pjj  -  2PSij/3)(l-C2)  __  ___ 

,J  P  +  pe(Ci-l)  f  Rik  =  *  2t0p(  EjpjpUku]  -  ekpjpu-u" ) 


3uj 


••  3uj 


Pij^-pu^-pUjU^) 


2P  =  Pii 


[2.14] 


where  Ci  =  1.5,  C2  =  0.6. 


In  the  near  wall  region  (viscous  sublayer  and  overlap  regions),  the  low-Reynolds 
number  form  employed  previously  is  used.  This  hybrid  approach  has  the  advantage  that  no 
near  wall  modelling  is  explicitly  required  for  the  ARSM  used,  but  the  influence  of  extra 
strain  rates  on  Reynolds  stress  distributions  are  accommodated. 

The  matching  point  for  the  hybrid  model  is  chosen  as  y^atch  =  200.  The 
justification  for  this  choice  is  as  follows.  As  pointed  out  by  Patel  et  al.  (1985),  the  near 
wall  shear  stress  damping  provided  by  the  functions  in  low  Reynolds  number  k-e 
formulations  are  designed  to  model  the  direct  influence  of  molecular  viscosity  on  the 
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Reynolds  shear  stresses.  However,  since  Reynolds  stresses  are  also  damped  due  to  the 
pressure-strain  correlation  (<J)‘j.w  above),  and  the  damping  constants  in  the  function  are 
obtained  by  empirical  matching  of  simple  shear  flow  data  (Chien  (1982)),  f^  implicitly 
accounts  for  near  wall  pressure  strain,  though  it  is  only  appropriately  correlated  for  the 
viscous  damping  effects,  in  light  of  this  and  the  choice  in  the  present  work  to  neglect  near 
wall  pressure  strain  modelling  in  the  ARSM  used,  y^match  was  chosen  as  200  since  for 
Chien's  model,  f^  has  nearly  reached  its  freestream  value  of  1.0  at  this  location 
(fji(y+  =  200)  =  .90). 


The  effective  stress  tensor  then  becomes  for  the  ARSM  : 


. ...  2  s  duk 

*■*  1  \9xj  3xJ  3  1J  3xfc 


-  pkTy  -  -2.  Sijpk 

3  [2.15] 


Equation  [2.15]  is  used  to  evaluate  viscous  fluxes  in  the  momentum  and  energy 
equations.  The  gradient  diffusion  approximation,  equation  [2.7],  is  retained  for  the  heat 
flux  terms  in  the  energy  equation.  Stability  implications  and  validation  of  the  hybrid  model 
are  provided  in  Chapters  3  and  4,  respectively. 

2.1.5  Transformation  of  Governing  Equations  to  Generalized  Coordinates 

In  general,  three-dimensional  flow  problems  of  engineering  interest  have  physical 
boundaries  which  can  be  very  complex.  It  is  convenient  in  these  circumstances  to  map  the 
physical  coordinates  to  a  rectilinear  computational  space,  (x,  y,  z)  -->  (£,  r|,  Q,  as  shown 
in  Figure  2.1.  In  this  thesis,  such  a  transformation  to  generalized  body  fitted  coordinates  is 
applied  to  the  seven  governing  equations  [2.3],  [2.10].  The  resulting  coupled,  Favre 


averaged,  three-dimensional,  continuity,  momentum,  energy  and  k  and  e  equations  can  be 
written  in  conservative  form  in  generalized  body  fitted  coordinates  as  : 
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Components  of  the  effective  stress  tensor,  tjj,  are  defined  in  equations  [2.8]  and  [2.15]. 
Metric  terms  in  equations  [2.17]  are  defined  by 
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In  this  thesis,  equation  [2.16],  is  solved  numerically  subject  to  appropriate  boundary 
conditions,  as  detailed  in  the  next  chapter.  For  notational  simplicity,  averaging  superscripts 
will  be  dropped  hereafter. 
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CHAPTER  3 

DEVELOPMENT  AND  ANALYSIS  OF  NUMERICAL  PROCEDURES 

3.1  Discretized  Governing  Equations 

The  numerical  solution  of  the  system  of  equations  [2.16]  subject  to  appropriate 
boundary  conditions,  is  the  topic  of  this  chapter.  The  spatial  and  temporal  discretization 
procedures  adopted  are  presented  first  An  extensive  discussion  of  artificial  dissipation, 
including  accuracy  and  stability  implications  on  highly  stretched  grids,  and  the  scalings 
incorporated  is  provided  in  section  3.2.  Boundary  and  initial  conditions  for  the 
turbomachinery  flow  computations  of  interest  in  this  thesis  are  presented. 

In  the  third  section  of  this  chapter,  a  comprehensive  numerical  stability  analysis  of 
the  discrete,  coupled  system  of  seven  governing  equations  is  presented.  Order  of 
magnitude  arguments  are  presented  for  flow  and  geometric  properties  typical  of  internal 
flows,  including  turbomachinery  applications,  to  ascertain  the  relative  importance  of  grid 
stretching,  rotation  and  turbulence  source  terms  and  effective  diffusivity  on  the  stability  of 
the  scheme.  It  is  demonstrated  through  both  analysis  and  corroborative  numerical 
experiments  that  1)  it  is  quite  feasible  to  incorporate,  efficiently,  a  two  equation  k-e 
turbulence  model  in  an  explicit  time  marching  scheme,  provided  certain  numerical  stability 
constraints  are  enforced;  2)  the  role  of  source  terms  due  to  rotation  on  the  stability  of  the 
numerical  scheme  is  negligible  for  rotation  numbers  typical  of  turbomachinery  flow 
calculations;  3)  the  direct  role  of  source  terms  in  the  turbulence  transport  equations  on  the 
stability  of  the  explicit  numerical  scheme  is  negligible,  except  in  the  earliest  stages  of 
iteration  (a  result  which  is  contrary  to  that  generally  perceived);  4)  there  is  no  advantage  to 
numerically  coupling  the  two  equation  model  system  to  the  mean  flow  equation  system,  in 
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regard  to  convergence  or  accuracy;  5)  for  some  flow  configurations,  including 
turbomachinery  blade  rows,  it  is  useful  to  incorporate  the  influence  of  artificial  dissipation 
in  the  prescription  of  a  local  timestep;  and  6)  explicit  implementation  of  an  algebraic 
Reynolds  stress  model  (ARSM)  is  intrinsically  stable  provided  that  the  discrete  two- 
equation  transport  model  which  provides  the  necessary  values  of  k  and  e  is  itself  stable. 


3.1.1  Spatial  Discretization.  Flux  Evaluation 


A  finite  difference  discretization  (for  evaluation  of  inviscid  and  viscous  fluxes), 
which  avoids  the  metric  singularity  at  the  interface  between  periodic  and  solid  wall 
boundaries  (i.e.  H-grid  leading  and  trailing  edges)  that  occurs  in  standard  finite  difference 
approaches,  has  been  adopted.  This  scheme,  presented  here,  removes  the  ambiguity  in 
metric  evaluation  at  such  boundaries,  and  results  in  a  more  compact  differencing  molecule 
at  all  grid  points.  In  the  present  discussion,  a  two-dimensional  formulation  is  presented  for 
the  purpose  of  brevity  and  clarity. 


Consider  a  two-dimensional  flux  evaluation  in  the  computational  plane  shown  in 

3Ec  3Ev 

Figure  3.1.  When  evaluating  fluxes  in  the  %  direction,  the  derivatives  ,  may  be 

discretized  to  second  order  as 


3E 


—  =  Ei+i/2,j  -  Ej.i/2,j 
or  ^ 


[3.1] 


where  ^  -  1.  The  former  gives  rise  to  a  13  point  differencing  molecule  in  two- 
dimensions  (33  point  in  three-dimensions)  for  the  viscous  fluxes. 
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Figure  3.1  Two-dimensional  representation  of  periodic  leading  edge  grid 
topology  in  physical  and  computational  space.  Fluxes  evaluated  at  X.  Dashed 
lines  represent  grid  in  periodic  adjacent  passage. 
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Referring  again  to  Figure  3.1,  evaluation  of  >+i/2.j  and  >+1/2J  require  values 
of  the  transport  variables,  pressure  and  effective  diff  usiv ities  at  i+1/2.  These  are  obtained 
unambiguously  by  averaging  values  at  i  and  i+1.  However,  if  the  metric  terms  required  at 
i+1/2  are  averaged  in  the  same  fashion,  the  metric  ambiguity  problem  mentioned  above 
manifests  itself.  To  see  this  consider  the  evaluation  of  ^  i+i/2.j  at  ile-1  as  illustrated  in 
Figure  3.1.  Straightforward  averaging  of  metrics  requires  values  of  £x,  T)x  and  T|y  at 
ile,  where 


=  Jyn  .  =  -Jxt,  ,  t\x  =  -Jy^  ,  riy  =  Jx^  ,  J  = 


_L 


xqyn  -  xny^ 


[3.2] 


Here,  x^  and  y^  are  not  uniquely  defined  at  i  =  ile,  because  either  the  suction  or  pressure 
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side  branch  can  be  considered  in  evaluating  ^ .  This  difficulty  is  avoided,  in  the  present 
technique,  in  that  all  metrics  used  in  flux  computations  are  evaluated  halfway  between  grid 
points  in  the  computational  domain.  Once  again  referring  to  Figure  3.1,  at  i+1/2,  we  have 
the  metric  quantities  as  aefined  in  equation  [3.2],  where  discretization  at  i+1/2  is  given  as  : 

i+1/2, j  =  xi+l,j  *  xi,j  j  y^i+l/2,j  =  yi+l.j  '  yi,j 
XT]  i+l/2,j  =  .25*(Xi+itj+i  *  xi+l,j-l  +  xi,j+l  '  xi,j-l) 
y-n  i+i/2j  =  .25*(yi+i,j+i  *  yi+i,j-i  +  yi.j+i  *  yij-i)  [3.3] 

Expressions  similar  to  equation  [3.3]  are  used  to  evaluate  fluxes  in  the  t\  and  £  directions  at 
j+1/2  and  k+1/2  respectively.  This  discretization  is  second  order  accurate.  Extension  of 
this  procedure  to  three  dimensions  is  straightforward. 

The  present  approach  requires  more  storage  for  metric  terms  than  if  the  fluxes  are 
computed  on  the  grid  vertices  themselves  (13  point  molecule),  however,  the  accuracy. 


convenient  application  of  boundary  conditions  and  numerical  stability  of  the  scheme,  all 
suffer  by  using  the  larger  differencing  molecule  when  evaluating  dissipation  fluxes,  so  only 
the  second  form  in  equation  [3.1]  was  considered.  Specifically,  the  truncation  error 
associated  with  discretizing  the  viscous  fluxes  on  a  uniform  cartesian  grid  using  a  thirteen 
point  scheme  is  0(4Ax2)  +  0(4 Ay2)  as  compared  to  0(Ax2>  +  0(Ay2>  for  the  nine  point 
molecule,  consistent  with  the  truncation  error  of  the  convective  fluxes.  It  is  also  easier  and 
more  consistent  to  apply  cyclic  differencing  along  periodic  boundaries  and  Neumann 
conditions  near  the  intersection  of  periodic  and  solid  boundaries  with  the  more  compact 
differencing  molecule.  In  addition,  it  can  be  shown  (Merkle  (1988))  that  using  the  larger 
differencing  molecule  is  equivalent  to  using  the  more  compact  form  with  added  negative 
dissipation,  that  is,  a  destabilizing  term. 

3.1.2  Explicit  Numerical  Procedure 

The  numerical  technique  used  for  temporal  discretization  in  the  present  studies 
incorporates  a  standard  4-stage  Runge-Kutta  scheme  as  first  applied  to  Euler  calculations 
by  Jameson  et  al.  (1981), 


Q  =  Q  +  ^  At  R(Q  ) 

""sri  1  i 

Q  =Q  +  j  At  R(Q  ) 

Q  =Q  +  ^-  At  R(Q  ) 

~n+l  C;n  .  ;>  .^3. 

Q  =  Q  +  At  R(Q  )  [3.4] 


Here,  the  residual,  R,  is  defined  according  to 
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Viscous  and  source  terms  are  evaluated  prior  to  the  first  stage,  convective  terms  are 
computed  at  every  stage.  The  scheme  is  first  order  accurate  in  time.  The  stability  region 
for  the  scheme  is  shown  in  Figure  3.2. 

The  curve  in  Figure  3.2  represents  the  contour  lg(z)l  =  1.0,  where  z  is  the  value  of 
the  complex  Fourier  symbol  of  any  discretized  scalar  equation  at  a  particular  wave  number, 
and  g  is  the  amplification  factor  arising  from  a  1-D  scalar  VonNeumann  linear  stability 
analysis  of  the  given  scheme  applied  to  this  discretized  equation. 

In  Figure  3.3,  an  amplification  factor  vs.  error  phase  plot  is  presented  for  this 
scheme  applied  to  a  scalar  advection  equation  with  a  CFL  number  of  2fl.  Also,  for 
comparison,  the  g  vs.  <J>  plot  for  a  Lax-Wendroff  scheme  (1960)  is  presented  where  a 
typical  CFL  number  of  0.85  was  used.  These  CFL  numbers  are  representative  of  those 
used  in  three-dimensional  Navier-Stokes  computations  using  the  Lax-Wendroff  scheme 
(Davis  (1991) )  and  Runge-Kutta  schemes.  This  diagram  illustrates  that  the  Runge-Kutta 
scheme  provides  no  high  wave  number  damping,  whereas  the  Lax-Wendroff  scheme 
intrinsically  damps  high  wave  number  errors.  This  illustrates  one  of  the  characteristics  of 
Runge-Kutta  procedures.  Specifically,  artificial  dissipation  must  be  purposefully  added  to 
the  scheme  for  stability,  conceptually  giving  the  user  more  control  of  its  levels. 

To  accelerate  the  solution  to  steady  state,  locally  varying  timesteps  are  computed 
based  on  a  linear  stability  analysis  of  the  discretized  Navier-Stokes  equations.  This  linear 
stability  analysis  of  the  discrete  coupled  7x7  system,  is  presented  in  section  3.3  of  this 
chapter,  and  provides  the  following  local  timestep  specification  : 
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Figure  3.2  Stability  boundary  of  Runge-Kutta  scheme  used  in  this  thesis.  The  curve 
represents  the  contour  lg(Z)l  =  1.0,  where  Z  is  the  complex  Fourier  symbol,  Z  =  X  + 
iY  of  any  discretized  scalar  equation  at  a  particular  wave  number,  and  g  is  the 
amplification  factor  arising  from  a  1-D  scalar  VonNeumann  linear  stability  analysis  of 
the  given  scheme  applied  to  this  discretized  equation. 


Figure  3.3  Amplification  factor,  Igl,  vs.  error  phase,  0,  plot  for  RK4  scheme  applied 
to  scalar  advection  equation.  For  comparison,  the  Igl  vs.  0  curve  for  the  Lax-Wendroff 

scheme  is  also  included 
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At  =  MIN  [  Ate  ,  Atv  ]  = 


-  16(c^  +  Gn  +  a;)  ] 


Here,  CFLop  and  Qop  arc  input  parameters  corresponding  to  operational  CFL  and 
VonNeumann  numbers  chosen  to  ensure  stability,  and  o^,  and  arc  general  artificial 
dissipation  scale  factors,  several  possible  forms  of  which  are  discussed  in  the  next  section. 
The  first  term  in  the  brackets  in  equation  [3.6]  arises  from  the  convection  operators,  the 
latter  terms  correspond  to  physical  and  artificial  dissipation.  The  influence  of  machine 
rotation,  turbulence  quantities  and  nonorthogonality  metrics  do  not  appear  in  this  stability 
specification,  as  justified  by  an  order  of  magnitude  analysis  carried  out  in  section  3.3  of 
this  chapter. 


Along  blade  surfaces,  the  no-slip  condition  is  imposed  upon  the  velocities, 
pressure  is  extrapolated  from  adjacent  grid  points,  and  density  is  computed  based  on 
specified  wall  temperature  or  heat  transfer  rate.  Also,  turbulent  kinetic  energy  and 
dissipation  rate  are  set  to  zero  along  solid  boundaries,  as  discussed  in  section  2.1.3. 
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Considering  the  %  direction  only,  the  inviscid  form  of  equation  [2.16],  can  be  written,  in 

the  absolute  frame,  in  nonconservative  form  as  ^  ,  where  A  is  the  7x7  flux 

-  A  =  ^ 

Jacobian  of  vector  be,  oy .  The  eigenvalues  of  A  are 

U,  U,  U,  U,  U,  U  +  aVv^- V£,  U  -  aVv^-  This  requires  that  for  subsonic  inlet 
flow  sue  boundary  conditions  should  be  specified.  Likewise,  for  subsonic  exit  flow  one 
boundary  condition  should  be  specified. 


Accordingly,  at  subsonic  inflow  boundaries,  total  temperature  profile  and  constant 

values  of  turbulent  kinetic  energy  and  dissipation  rate  are  specified.  Following  Kirtley  et 

al.  (1990),  for  the  subsonic  rotor  flow  computation  presented  in  Chapter  5,  and  for  the 

transonic  centrifugal  compressor  computation  presented  in  Chapter  6,  inlet  flow  angles  and 

the  axial  component  of  mass  flux  are  specified,  and  the  static  pressure  is  extrapolated  along 

T)  =  constant  grid  lines  from  the  interior  of  the  computational  domain.  This  allows  the  user 

to  impose  a  desired  operating  flow  coefficient  almost  exactly,  without  having  to  adjust  exit 

static  pressure  so  as  to  pass  the  appropriate  mass  flow  to  match  the  flow  coefficient.  For 

all  other  calculations,  either  flow  angle  or  crossflow  velocity  components  and  total  pressure 

R-  =  V  -  -2e- 

profile  are  specified,  and  the  R'  characteristic,  Y  *  * ,  is  extrapolated  along  T|  = 

constant  grid  lines  from  the  interior  of  the  computational  domain.  Specifically,  designating 
subscripts  1  and  2  for  grid  planes  at  and  adjacent  to  the  inlet  boundary : 


_(Y-  1)(R'2)2  +  V2(1-Y)(R2)2  +  4(y+  l)cpTo2 
(y+  1) 


a  result  which  is  easily  obtained  using  the  definitions  of  total  temperature  and  the  left 
running  Riemann  variable. 
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At  subsonic  outflow  boundaries,  static  pressure  is  specified  and  relative  velocity 
components,  entropy,  pk  and  pe  are  extrapolated  along  r|  =  constant  grid  lines.  For 
annular  geometries,  a  simplified  radial  equilibrium  equation  : 

dp  =  pve2 

dr  P  r  [3.8] 

is  used  to  integrate  a  specified  hub  exit  static  pressure  in  the  radial  direction  to  provide  the 
exit  pressure  distribution.  Along  periodic  boundaries,  values  of  the  flowfield  variables  arc 
prescribed  along  periodic  "ghost"  grid  points  and  this  information  is  used  when  discretizing 
derivatives  in  the  T|  direction. 

The  constant  values  of  pk  and  pe,  which  are  imposed  at  the  inflow  boundary,  are 
based  on  specified  freestream  turbulence  intensity  and  length  scale. 


Typically ,  the  freestream  length  scale  is  set  between  .001  and  .01  times  the  blade 
spacing,  as  this  scale  is  of  the  same  order  as  the  largest  eddy  motions,  i.e.  the  thickness  of 
the  boundary  layers  which  develop  along  the  blades  and  endwalls. 

The  flowfield  is  initialized  using  an  inviscid  quasi- ID  analysis  to  provide  uniform 
initial  velocity  profiles  along  each  ^  =  constant  grid  line.  In  situations  where  this  simple 
initialization  strategy  has  proven  to  be  insufficient,  alternative  procedures  are  detailed  in  the 
following  chapters. 


When  central  difference  schemes,  which  do  not  inherently  damp  high  wave 
number  errors,  are  applied  to  hyperbolic  problems,  the  addition  of  artificial  dissipation  is 
required  to  damp  high  wave  number  disturbances.  These  disturbances  can  be  introduced 
into  linear  problems  through  inconsistent  boundary  condition  treatment  or  machine 
roundoff  error.  In  non-linear  problems,  such  disturbances  can  be  introduced  through 
aliasing  of  sub-grid  scale  non-linear  disturbances,  to  lower,  resolvable  wave  numbers. 
Even  for  viscous  flow  calculations,  artificial  dissipation  must  be  introduced  into  the 
scheme  because  the  physical  viscous  terms  are  only  effective  in  damping  frequencies  at 
higher  wave  numbers  than  can  be  resolved  on  practical  grids. 


In  the  present  work,  artificial  dissipation  is  added  to  the  five  discretized  mean  flow 
equations  [3.5]  as : 


/  ^ 
3Ec 

,3FC 

3GC\ 

U 

<*1 

as) 

D(Q)  +  S 


[3.10] 


Here,  D(Q)  represents  a  mixed  2nd  and  4th  order  nonconservative  artificial 
dissipation  operator  similar  to  that  devised  by  Jameson  et  al.  (1981) : 

D(Q)  =  D^(Q)  +  DT1(Q)  +  D;(Q) 

D^(Q)  =  S2^Q  +  S^5^Q 


D-q(Q)  —  S2t|$t|tiQ  S4T)§T|Trr|T)Q 
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D;(Q)  =  S2;5^Q  +  [3  j  j 

The  fourth  order  operators  are  included  to  damp  high  wave  number  errors  and  the 
second  order  operators  are  included  to  improve  shock  capturing. 

As  pointed  out  by  Pulliam  (1986),  artificial  dissipation  terms  should  operate  on 
physical  values  of  the  flowfield  variables,  that  is,  D(Q),  as  in  equations  [3.1 1],  not  D(Q). 
This  is  required  since  in  general,  regions  of  smoothly  varying  flowfield  properties  can 
occur  in  regions  of  non-smooth  grid  distribution.  Accordingly,  the  scale  factors,  S,  in 
equations  [3.1 1]  should  be  scaled  by  the  inverse  of  the  metric  Jacobian,  J. 

In  addition  to  the  above  consistency  requirement  on  the  dissipation  scaling, 
artificial  dissipation  should  always  be  reduced  to  levels  adequate  to  stabilize  a  scheme 
without  altering  the  accuracy  of  the  solution.  For  the  computation  of  viscous  flows  on 
highly  stretched  grids,  this  requirement  is  a  sensitive  matter. 


If  flux  split  differencing  is  used  for  the  convection  operators,  an  inherent  amount 
of  artificial  dissipation  is  introduced.  It  can  be  easily  shown  that  flux  split  differencing 
schemes  are  equivalent  to  central  differencing  schemes  with  appropriately  added  artificial 
dissipation.  Since  Runge-Kutta  schemes  are  intrinsically  less  dissipative  than  upwind 
schemes,  it  is  attractive  to  purposefully  add  artificial  dissipation  to  these  schemes  in  a  way 
as  to  provide  the  stability  and  shock  capturing  advantages  of  upwind  schemes,  while 
retaining  control  of  its  levels.  Accordingly,  flux  split  approaches  can  be  used  to  "guide" 
the  development  of  central  difference  artificial  dissipation  schemes  [Pulliam  (1986)], 
leading  to  what  is  referred  to  in  this  thesis  as  eigenvalue  scaling. 
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By  way  of  example,  Pulliam  (1986)  shows  that  second  order  Steger-Warming 

3  Eg 

(1981)  flux  vector  splitting  applied  to  the  convection  operator ,  ,  is  equivalent  to  a 

central  difference  plus  a  fourth  order  difference  operating  on  Q.  The  coefficient  of  this 
operator  is  the  flux  Jacobian  matrix,  A  (where  Ec  =  AQ).  This  suggests  the  use  of  the 
spectral  radius  of  matrix  A  as  a  scale  factor  in  artificial  dissipation  operators  applied  to 
central  difference  schemes. 

In  explicit  local  time  stepping  schemes  applied  to  inviscid  problems,  local 
timesteps  are  generally  computed  such  that  they  are  inversely  proportional  to  the  spectral 
radius  of  A,  p(A),  or  in  multi-dimensions,  some  linear  combination  of  the  spectral  radii  of 
the  several  flux  Jacobians.  For  this  reason,  most  early  workers  (see  Jameson  et  al.  (1981) 
for  example)  scaled  dissipation  operators  by  p(A)  +  p(B)  (for  two-dimensional  Euler 

±  At  =  --lA-ty- 

computations;  Fc  =  BQ),  or  equivalently,  by  At  where  +  Aty .  Here,  Atx  and  Aty 

are  inversely  proportional  to  p(A)  and  p(B)  respectively.  From  a  flux  split  point  of  view, 
this  says  that  the  artificial  dissipation  is  proportional  to  the  propagation  speed  of  the 
characteristic  which  limits  the  maximum  stable  local  timestep.  This  approach  has  been 
widely  used  in  multidimensional  central  difference  inviscid  flow  computations  where  grids 
are  relatively  uniform. 

For  high  Reynolds  number  flows,  very  highly  stretched  grids  must  be  used  to 
resolve  body  normal  gradients  in  near  wall  regions.  If  the  artificial  dissipation  terms  in 
both  the  ^  and  tj  directions  are  scaled  by  the  local  timestep  on  grids  which  are  highly 
stretched  in  the  ri  direction,  excessive  dissipation  is  introduced  in  the  £  direction.  This 
effect  is  discussed  by  Caughey  and  Turkel  (1988)  and  Swanson  and  Turkel  (1987).  This 
excessive  dissipation  may  reduce  accuracy  and  convergence  rates  in  viscous  flow 
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computations.  Martinelli  (1987)  recently  devised  a  flux  Jacobian  eigenvalue  scaling  of  the 
artificial  dissipation  terms  which  alleviates  this  problem.  Since  the  present  technique  is 
primarily  used  to  compute  viscous  flows  on  highly  stretched  grids,  anisotropic  dissipation 
scaling  factors  similar  to  those  used  by  Martinelli  are  incorporated.  In  his  two-dimensional 
work,  Martinelli  provided  the  following  form  for  the  scale  factors  in  equation  [3. 1 1] 

(where  the  present  finite  difference  notation  is  used) : 


S25  = 


s4^  = 
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[3.12] 


Here,  Atc^,  At^,  are  timesteps  corresponding  to  unit  CFL  limit  for  the  inviscid 
one-dimensional  problem  in  each  direction. 


Atd;  = 


|U|  +  aV(v^V$) 


Atcr|  — 


|V]  +  aV[V"rvVrJ 


[3.13] 


and  o^,  02^,  o4^,  cj4t1  are  scalar  coefficients  discussed  below. 

The  choice  of  unit  CFL  scaling  in  the  numerator  of  equation  [3.13]  ensures  that  the 
steady  state  solution  will  be  independent  of  the  operational  CFL  number  used  to  compute 
local  timesteps.  If  a  =  1,  equation  [3.12]  reduces  to  standard  isotropic  scaling.  As 
mentioned  above,  this  introduces  excessive  dissipation  in  the  4  direction  in  regions  where 
the  grid  is  stretched  in  the  T)  direction.  If  a  =  0,  the  scaling  becomes  purely  anisotropic.  If 
the  grid  is  very  highly  stretched  in  the  T)  direction,  such  scaling  may  not  provide  enough 
dissipation  in  the  %  direction,  resulting  in  reduced  convergence  rates.  For  intermediate 
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values  of  a  between  1/2  and  2/3,  Martinelli  (1987)  and  Swanson  and  Turkel  (1987)  have 
shown  good  convergence  rates  for  Euler  and  Navier-Stokes  calculations  on  highly 
clustered  grids. 

Careful  examination  of  the  scaling  in  equation  [3.12],  shows  that  if  the  coefficient, 
a,  is  reduced  to  zero  (towards  purely  anisotropic  scaling),  local  levels  of  artificial 
dissipation  do  decrease  in  the  direction  opposite  that  in  which  the  grid  is  clustered, 
consistent  with  the  arguments  presented  above.  However,  the  dissipation  levels  in  the 
direction  of  clustering  actually  increase.  To  see  this,  consider  a  two-dimensional  grid 
element,  stretched  in  the  q  direction,  of  aspect  ratio  100.  If  this  element  lies  within  the 
boundary  layer,  where  local  velocity  is  small,  the  spectral  radii  of  the  flux  Jacobians  can  be 
approximated  by : 

p(A)  »  |U|  +  aVv^-V^  =  ar/v\- VE,,  p(B)  =  |V[  +  aVVq-Vq  =  aVVq-Vq 
The  scaling  [3.12]  provides  for  a  =  1  (isotropic) : 

=  aVv^-V^  +  aVVq-Vq  =  100  K 
Sn  =  aVv^- +  aVVq- Vq  =  100  K 

where  K  is  some  constant,  proportional  to  a  constant  length  scale  in  the  physical  domain 
and  the  local  speed  of  sound.  If  a  purely  anisotropic  scaling  is  introduced,  a  =  0  : 

S;:  =  2aV  V  ^  •  V  ^  =  2K 


St,  =  2aVVq-Vq  =  200  K 
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that  is,  doubles  in  magnitude,  which  substantiates  the  above  argument.  It  has  been  the 
experience  of  the  author,  that  in  three-dimensional  viscous  computations,  it  is  crucial  to 
incorporate  very  low  values  of  a,  in  order  to  obtain  stable  solutions  (  a  =  0. 1  was  used  for 
the  rotor  flow  computation  presented  in  Chapter  5,  for  example).  Accordingly,  Martinelli’s 
scaling,  equation  [3.12],  has  been  modified  and  extended  to  three-dimensions  in  this  work 
as  follows : 
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[3.14] 


where 


At<*  = 


|Uj  +  aVJvVV^J 


,  Ateq  — 


|V|  +  aV{Vr|- VrJ 


-,  Atcr  - 


[W|  +  aV(v;-VC’ 


If  this  scaling  is  applied  to  the  example  provided  above,  for  a  =  1  : 


[3.15] 


=  aVv^- +  aVv-q- Vr)  =  100  K 
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=  aVv^- VE,  +  aVVr|- Vrj  =  100  K 
identical  to  the  original  scaling,  but  for  a  =  0  : 


=  aV  V  ^  •  V  E,  =  K 


=  aVVr]- Vri  s  100  K 

That  is,  the  modified  anisotropic  eigenvalue  scaling,  does  not  increase  levels  of 
dissipation  in  the  direction  of  grid  stretching.  This  modified  scaling,  published  originally 
in  Kunz  and  Lakshminarayana  (1990)  has  been  adopted  by  Chima  (1991),  who  has  also 
found  this  scaling  to  provide  a  "good  distribution  of  dissipative  terms  in  each  direction”,  in 
three-dimensional  viscous  flow  calculations. 

If  no  shocks  are  anticipated  in  the  flowfield,  the  scalar  coefficients,  a,  in  equations 
[3.14]  are  defined  as  o4^  =  o4ti  =  o4^  =  -  k4,  where  typically  =  0.02.  If  shocks  are 
anticipated  in  the  flowfield  (see  the  supersonic  cascade  results  in  Chapter  4),  a  composite 
dissipation  operator  devised  by  Jameson  et  al.  (1981)  is  used  to  turn  off  locally  the  fourth 
order  central  difference  dissipation  and  turn  on  a  second  order  central  difference  dissipation 
operator.  In  this  case,  the  scalar  coefficients,  a.  in  equations  [3.14]  take  the  following 
form : 


-  K2  MAX  (Vr1  +  li  j,  k,V,,.  J<  k.V^M,  J,  k) 


04^  —  -  max  (i»,  -  o 2£ j 


13.16] 


where  the  dissipation  "sensor"  v,  is  a  normalized  second  derivative  of  pressure  that  is  very 
small  except  in  the  immediate  vicinity  of  shocks  and  stagnation  points  : 


v  .  .  _  |Ph-1,  j,  k  -  2  Pi.  j.  k  +  Pi-1,  j,  k) 

^l,J'  Pi+1,  j.  k  +  2  Pi,  j,  Ic  +  Pi-1,  j.  k  [3.17] 

Here,  k2  =  0.25,  k4  =  0.02.  Expressions  similar  to  equations  [3.16]  and  [3.17]  are  used 
in  the  T|  and  £  directions.  When  shocks  are  not  anticipated  in  the  flowfield,  k2  is  set  equal 
to  zero  so  that  fourth  difference  artificial  dissipation  only  is  added  to  the  mean  flow. 


The  dissipation  scaling  issues  presented  so  far  are  essentially  inviscid  in  nature,  in 
that  eigenvalues  of  the  inviscid  flux  Jacobians  only  are  considered.  As  the  computational 
grid  in  a  Navier-Stokes  solution  is  refined,  smaller  and  smaller  scales  can  be  resolved. 
Eventually,  the  small  length  scales  associated  with  nonlinear  disturbances  introduced  by 
convection  operators  in  the  governing  equations  (shocks  for  instance)  can  be  resolved,  and 
the  gradients  that  these  disturbances  give  rise  to  can  be  effectively  damped  by  physical 
dissipation.  Accordingly,  it  is  desirable  that  any  variable  artificial  dissipation  scaling 
should  have  the  property  that  as  the  grid  is  refined,  artificial  dissipation  goes  to  zero.  This 
item  is  considered  here;  again  a  two-dimensional  formulation  is  adopted  for  clarity. 


2 

If  a  cartesian  physical  grid  is  considered,  with  Ax  =  Ay,  ^ =  ^ , 

VtV  Vrl  =  "Hy  (£  coincident  with  x-axis,  T)  coincident  with  y-axis),  then  the  scaling  factors 
proposed  in  equations  [3.14]  take  on  the  following  forms  for  isotropic  (a  =  1,  subscript  i) 
and  anisotropic  (a  =  0,  subscript  a)  scalings  : 


|Uj  +  |V]  +  aV  V  ^  •  V  ^  +  aVVq-  Vq 


=  Ayo(|iij+a)  +  Axo(M+a) 


S?t,  =  Ayo(|uj+a)  +  Axoflvj+a) 


=  Ayo(|ul+a) 
S2r|  =  Ax0(H+a) 


[3.18] 


where  superscript  0  denotes  the  grid  under  consideration.  If  the  grid  is  doubled  n  times  in 
the  rj  direction,  with  n  »  1,  as  would  be  representative  near  a  solid  boundary  in  a  viscous 
flow  computation,  then  denoting  the  refined  grid  by  superscript  1,  the  scale  factors  become 

Sfc  =  AyiOu|+a)  +  Axo{|vj+a)  =  ^-Ayo(|u|+a)  +  Axoflv|+a)  =  Axoflvj+a) 

S.'ti  2  Axoflv|+a} 

=  ^-Ayoflul+a) 

Sin  =  Axo(jv|+a)  [3.19] 

So  clearly,  as  the  grid  is  refined  ;  i  one  direction  only,  both  isotropic  and 
anisotropic  scale  factors  remain  finite  as  the  grid  is  refined  infinitely.  However,  " 

>  0  as 

Ay  — >  0,  because,  unlike  the  fourth  derivative  of  a  function  Q,  which  may  be  large  if  Q  is 
exponential  in  nature,  the  fourth  difference  stencil  applied  to  Q  will  always  tend  to  zero  as 
the  grid  is  refined  for  any  continuous  function  Q.  So  the  eigenvalue  scaled  artificial 
dissipation  operators  do  have  the  desirable  property  that  artificial  dissipation  goes  to  zero  as 
the  grid  is  refined,  even  if  in  only  one  direction,  as  is  typical  in  high  Reynolds  number 
viscous  flow  computations. 

A  problem  arises  in  viscous  flow  computations  though,  because  the  rate  at  which 
artificial  dissipation  goes  to  zero  with  grid  refinement  may  not  be  adequate  to  keep  levels  of 
artificial  dissipation  significantly  less  than  the  physical  dissipation  very  near  a  solid 
boundary.  This  is  probably  not  surprising  since  all  of  the  scalings  presented  are  based  on 
inviscid  considerations  only.  This  observation  warrants  additional  boundary  layer  scaling 


such  that  artificial  dissipation  remain  significantly  less  than  physical  dissipation  through 
most  of  the  boundary  layer. 
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Another  scaling  issue  is  important  in  the  computation  of  viscous  flows.  Near  solid 
boundaries,  the  viscous  fluxes  in  the  momentum  equation  are  quite  large  and  are 
themselves  adequate  to  provide  smoothing.  As  discussed  above,  in  these  same  regions, 
second  and  fourth  derivatives  of  the  transport  variables  can  be  quite  large  leading  to  large 
values  of  artificial  dissipation  there.  This  phenomenon,  which  has  been  observed  by 
several  other  researchers  (Davis  et  al.  (1986)  and  Swanson  and  Turkel  (1987)  for 


instance),  can  give  rise  to  very  large  nonphysical  values  of  total  dissipation  in  the  near  wall 
region.  It  has  been  the  experience  of  the  present  author,  that  in  the  momentum  equations, 

/V  /N  ^  /S 

the  ratio  of  artificial  to  physical  dissipation,  i.e.  Rart/Rphys(  where  Ran  =  D(Q),  Rphys 

/aiv  3FV  J  3GV\ 

I - +  ~^rf  + - 

=1  ^  can  be  as  high  as  O(101),  in  the  very  near  wall  region  if  the  grid  is 

highly  stretched  locally!  The  accuracy  implications  of  this  are  self-evident.  Personal 
experience,  and  that  of  others  (Kirtley  (1990),  Weiss  (1991)),  have  shown  non-physical 
overshoots  of  nearwall  velocity  profiles  in  converged  viscous  flow  solutions  on  highly 
stretched  grids,  when  artificial  dissipation  operators  which  incorporate  inviscid  flow 
physics  in  their  formulation,  are  used. 


In  this  thesis,  a  geometric  decay  function  is  used  to  control  the  levels  of  artificial 
dissipation  in  near  wall  regions,  in  order  to  reduce  the  magnitude  of  numerical  to  physical 
smoothing  to  acceptable  levels.  Specifically,  the  scaling  functions  in  equations  [3.14]  are 
multiplied  by  a  normalized  square  of  the  local  velocity,  (W/Wrcf)2,  for  the  momentum 


equations  : 
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1  D(Q) 


[3.20] 


In  Figure  3.4,  this  scaling  function  is  plotted  as  a  function  of  W/Wref.  As  shown 
below,  this  additional  local  velocity  scaling  can  greatly  improve  solution  accuracy  for  the 
high  Reynolds  number  viscous  flow  computations  of  interest  in  this  thesis. 


To  illustrate  the  above  considerations  in  a  quantitative  fashion,  a  series  of 
numerical  experiments  were  undertaken. 


In  Chapter  4,  a  turbulent  flat  plate  flow  is  computed  to  validate  the  turbulence 
models  used.  The  details  of  these  computations  are  provided  there;  suffice  it  here  to  say 
that  the  Reynolds  number  based  on  distance  from  the  leading  edge  for  this  flow  reaches 
OflO6)  and  the  grid  is  therefore  very  refined.  If  a  typical  value  of  K4  =  0.02  is  used  with 
the  modified  eigenvalue  scalings  presented  in  equations  [3.14]  and  [3.16],  and  the  flat  plate 
flow  is  computed,  a  standard  law-of-the-wall  velocity  profile  similarity  is  not  achieved. 
This  is  illustrated  in  Figure  3.5a,  where  converged  boundary  layer  velocity  profiles  are 
plotted  in  inner  variables,  at  a  streamwise  location  where  the  local  Reynolds  number  based 
on  boundary  layer  thickness  is  approximately  80000  [to  match  Klebanoffs  (1954) 
experiments;  see  Chapter  4].  Clearly,  local  velocity  is  seen  to  undershoot  a  standard  log- 
law  profile. 

The  artificial  dissipation  scale  factor,  K4,  had  to  be  reduced  by  a  factor  of  40  in 
order  to  bring  the  converged  solution  to  within  acceptable  agreement  of  law-of-the-wall 
profiles.  This  is  illustrated  in  Figure  3.5b,  where  K4  has  been  reduced  to  0.0005. 


Figure  3.5  Effect  of  artificial  dissipation  velocity  scaling  on  law-of-the-wall  predictions  for  a 
turbulent  flat  plate  flow  computation.  Standard  law-of-the-wall  relations  (solid  lines)  and 
computed  values  (symbols).  Refer  to  section  4.1  for  details  of  configuration  and  computation. 
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If  free  stream  velocity  is  used  for  Wref,  and  equation  [3.20]  is  used  to  provide 
additional  scaling  in  near  wall  regions,  the  standard  value  of  K4  =  0.02  can  be  used  without 
corrupting  the  solution  accuracy.  This  is  illustrated  in  Figure  3.5c. 

It  is  noted  that  although  a  stable  solution  could  be  obtained  for  the  flat  plate  flow, 
using  K4  =  0.0005,  such  a  low  value  would  bring  about  a  rapidly  unstable  solution  in 
typical  engineering  calculations,  where  the  flow  physics  are  not  diffusion  dominated.  To 
illustrate  the  above  scaling  issue  for  such  a  typical  high  Reynolds  number  turbomachinery 
computation,  the  low  subsonic  cascade  flow,  described  in  detail  in  Chapter  4,  was  subject 
to  a  similar  numerical  experiment.  Figure  3.6  shows  the  converged  suction  surface  skin 
friction  coefficient  predictions  for  the  same  case  run  with  and  without  the  velocity  scaling 
defined  in  equation  [3.20].  This  diagram  again  illustrates  that  significant  errors  can  be 
introduced  due  to  near  wall  artificial  dissipation  levels  in  viscous  computations  on  highly 
stretched  grids. 

It  is  also  worth  noting  that  the  dissipation  scaling  considerations  addressed  above 
are  especially  important  when  a  multigrid  acceleration  scheme  is  used.  Careful  tuning  of 
artificial  dissipation  levels  is  crucial  when  performing  explicit  multigrid  calculations  on 
highly  stretched  grids.  This  is  because  inadequate  or  excessive  dissipation  can  diminish  the 
high  wave  number  damping  properties  of  the  driving  scheme  thereby  rendering  multigrid 
acceleration  less  effective  Swanson  (1990). 

Lack  of  adequate  grid  resolution  in  the  £  direction  just  upstream  and  downstream 
of  the  blade  edges,  causes  the  wall  damping  function,  f^,  in  equation  [2.12]  to  be  effective 

only  over  two  to  four  grid  points  in  this  direction.  This  gives  rise  to  very  large  streamwise 
gradients  in  k  and  e  at  leading  and  trailing  edges,  which  in  turn  leads  to  slowly  growing 
oscillations  in  these  variables.  It  has  therefore  been  found  necessary,  in  turbomachinery 


Figure  3.6  Comparison  of  predicted  skin  friction  coefficient,  along  the  suction  surface 
of  a  subsonic  double  circular  arc  cascade,  with  and  without  velocity  scaling  of  artificial 

dissipation. 
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blade  row  computations,  to  smooth  these  oscillations  by  incorporating  small  amounts  of 
second  order  artificial  dissipation  in  the  k  and  e  equations,  in  the  streamwise  direction  only, 
in  a  manner  consistent  with  equation  [3.14],  with  S2^  =  0.01,  S2ti  =  S2^  =  =  S4n  = 

S4?  =  0.0. 


Both  the  eigenvalue  scaling  (equation  [3.14])  and  local  velocity  scaling  (equation 
[3.20])  were  used  in  all  computations  that  follow.  The  precise  values  of  the  reference 
relative  velocity  used  in  the  scaling  equation  [3.20],  are  stated  for  each  of  the  computations 
presented. 


In  order  to  provide  some  verification  of  the  conclusions  of  the  analyses  provided  in 
this  section,  a  number  of  parametric  studies  have  been  undertaken  using  the  code  applied  to 
a  series  of  numerical  test  cases.  The  effects  of  rotation  and  turbulence  source  terms,  grid 
stretching  and  effective  diffusivity  and  artificial  dissipation  are  to  be  investigated,  and 
therefore  these  test  cases  should  embody  some  or  all  of  these  features.  Also,  the  test  cases 
should  typify  engineering  calculations  of  complex  internal  flowfields  such  as  occur  in 
turbomachines.  Accordingly,  a  three-dimensional  subsonic  compressor  rotor  and  a  two- 
dimensional  cascade  flow  configuration  are  used.  In  Chapters  4  and  5, 180,225  and 
1 1,739  point  grids  have  been  used  to  compute  these  flows,  where  the  method  is  seen  to 
capture  much  of  the  complex  viscous  flow  physics,  within  the  accuracy  of  the  turbulence 
and  tip  clearance  models  used.  The  computational  grids  and  configuration  details  for  these 
cases  are  included  there.  The  author  feels  that  these  grid  sizes  are  typical  of  those  used  for 
viscous  turbomachinery  flow  computations  in  industry.  Additionally,  for  some  of  the 
parametric  studies  performed  herein,  it  was  necessary  to  run  the  code  many  times.  Every 
fourth  grid  point  was  retained  from  the  "production"  3D  rotor  grid  for  these  cases.  The 
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resulting  grid,  is  much  too  coarse  to  accurately  resolve  the  near  wall  physics,  but 
convergence  properties  are  similar  to  those  of  the  fine  grid  solution. 

A  standard  four  stage  Runge-Kutta  scheme  was  chosen  for  its  familiarity  in  the 
numerical  verification  studies  that  follow,  but  the  results  obtained  are  valid  for  any  Runge- 
Kutta  scheme. 

3.3.1  Treatment  of  Mean  Flow  and  Turbulence  Transport  Systems 

Researchers  who  have  implemented  two-equation  turbulence  models  in  time¬ 
marching  Navier-Stokes  procedures,  have  often  adopted  schemes  where  the  turbulence 
transport  equations  are  numerically  decoupled  from  the  mean  flow  equations  (Grasso  and 
Speziale  (1989),  Holmes  and  Connell  (1989),  Liu  (1987),  Sahu  and  Danberg  (1986), 
Gerolymos  (1990),  Mavriplis  and  Martinelli  (1991)).  In  the  present  context,  numerical 
uncoupling  implies  that  at  each  iteration  the  mean  flow  equations  are  integrated  in  time, 
using  frozen  values  of  k  and  E  from  the  previous  iteration.  The  mean  flowfield  is  then 
frozen  and  the  turbulence  equations  are  integrated.  This  type  of  procedure  is  referred  to  as 
"uncoupled"  by  Grasso  and  Speziale  (1989),  but  elsewhere  in  the  literature  as  "lagged" 
(Sahu  and  Danberg  (1986))  or  "split"  (Lee  and  Dulikravich  (1991)).  In  explicit  multistage 
schemes,  "coupling"  simply  means  that  fluxes  are  evaluated  at  each  stage  using  values  of 
transport  variables  obtained  at  the  previous  stage  of  the  scheme  rather  than  using  some 
values  from  the  previous  iteration.  This  should  not  be  confused  with  the  standard  coupling 
of  systems  of  equations  in  many  implicit  schemes. 

Simple  linear  stabilty  analysis  indicates  that  the  stability  characteristics  of  the 
"coupled"  and  decoupled  approaches  are  identical.  In  earlier  work  by  the  present  author 
(Kunz  and  Lakshminarayana  (1990))  a  decoupled  approach  was  chosen  to  simplify  code 
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development,  and  was  used  for  the  two-dimensional  results  presented  in  Chapter  4.  Others 
have  cited  code  modularity  (Gerolymos  (1990))  and  the  ability  to  advance  the  equation  sets 
at  different  local  timesteps  (Holmes  and  Connell  (1989))  as  advantages  to  an  decoupled 
approach.  More  recendy,  a  coupled  approach  has  been  adopted  by  the  author;  all  three- 
dimensional  results  presented  in  Chapters  5  and  6  used  this  approach.  It  has  been  the 
experience  of  the  author  that  this  gives  rise  to  a  more  compact  and  better  organized  code 
and  has  no  affect  at  all  on  the  stability,  accuracy  or  CPU  performance  of  the  scheme. 

3.3.2  Vector  VonNeumann  Procedure 


In  generalized  curvilinear  coordinates,  the  compressible  three-dimensional 
continuity,  momentum,  energy,  turbulent  kinetic  energy  and  turbulent  kinetic  energy 
dissipation  rate  equations  can  be  written  in  conservation  form  as  : 


3Q  _  ai  3FC  38c  (3E,  .  8F,  .  3G,  \  a 

*‘1aT  ^  If  IT  *7 


[2.16] 


where  the  transport  variable,  flux  and  source  vectors  are  given  in  equations  [2.17].  When 

/N 

an  effective  diffusivity  formulation  is  adopted,  we  have  ~  ^(de.QjQ^QTi.Q;) , 

Fv  =  FvfpecQ.Q^Q^Q^) ,  Gv  =  Gv^cQ.Q^QtvQ;)  To  simpiify  the  problem,  viscous 
flux  vectors  are  linearized  as  ~  Fv(Q^Qn,Q<;)  Fv  -  FvfQ^.Q^Q^ 

Gv  -  GV(Q^,Q11,Q^)  ^  stability  analysis,  the  diffusivities  are  taken  as  locally 

^4  -SijPk) 

constant,  and  the  influence  of  the  components  of  the  normal  stress  terms,  ®x>  are 

ZP  .  VT 

neglected,  as  is  a  heat  transfer  term  prcpordonal  to  P  .  Such  linearizations  have 


also  been  incorporated  in  stability  analyses  performed  by  other  authors  (MacCormack  arid 
Baldwin  (1975),  Abarbanel  and  Gottlieb  (1981)).  These  assumptions  are  not  invoked  in 


These  flux  and  source  Jacobians  are  provided.  Matrix  P,  which  scales  the  viscous  and 
source  terms  consistent  with  the  operations  performed  to  convert  to  nonconservative  form, 
also  appears  in  Appendix  A.  In  the  linearization  of  source  Jacobian  D  ,we  take  S  =  S(Q). 
Such  a  linearization  has  be  .1  used  in  an  uncoupled  implicit  scheme  by  Sahu  and  Danberg 
(1986)  and  in  uncoupled  explicit  schemes  which  incorporate  pointwise  implicit  treatment  of 
the  turbulence  equation  source  terms  by  Liu  (1987)  and  Gerolymos  (1990).  Additionally, 
the  density  is  taken  as  locally  constant  in  linearizing  the  near  wall  modelling  terms  D  and  E. 
This  is  justified  since  these  functions  are  very  small  except  in  the  immediate  vicinity  of  the 
wall  where  local  Mach  numbers  are  low.  It  is  noted  that  viscous  deformation  work  terms 


do  not  appear  in  the  nonconservative  form  of  the  energy  equation  due  to  the  choice  of 
pressure  as  the  dependent  variable. 
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Fourth  order  difference  artificial  dissipation  operators  are  incorporated  in  the 
stability  analysis  (inclusion  of  second  order  dissipative  terms  follows  analogously  and  is 
not  included  for  brevity),  so  we  add  the  following  to  the  RHS  of  equation  [3.21]  : 


(equation  [3.21]) 
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[3.22] 


where  o^,  o^,  o^  are  any  general  locally  varying  dissipation  scale  factors.  These  factors 
are  real  and  negative  (so  that  the  fourth  difference  operators  be  dissipative)  and  have 
dimension  f1. 

If  second  order  accurate  central  differences  are  used  to  discretize  the  spatial 
derivatives,  the  corresponding  Fourier  symbols  for  the  spatial  derivatives  in  equations 
[3.21]  and  [3.22]  become  : 


Q^— »iS{: 

Q25-1 

Q^q-I  f 

Q2ti~ 

Q«r*-s$ss 

Q4t)~^4^Ct|-1)^ 

Q;-*iS  t; 
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where  ^  "  s’n^’  Q;  “  cos^\  etc...,  0  is  the  phase  of  the  error  mode,  and 
=  Aq  =  =  1  The  Fourjer  symbols  given  above  are  consistent  with  the  flux 

differencing  in  the  code.  This  is  an  important  consideration  as  factors  of  2  and  4  can 
appear  in  the  final  stability  expressions  if  viscous  fluxes  are  discretized  differently. 


If  a  standard  multi  stage  Runge-Kutta  scheme  is  incorporated,  then  a  vector 
VonNeuman  stability  analysis  yields  a  polynomial  expression  for  amplification  matrix,  G  = 
P(AtZ),  where  P  represents  some  matrix  polynomial  function,  and  Z  is  the  Fourier  matrix 


66 


of  the  spatial  discretization  of  the  governing  equations.  A  necessary  condition  for  stability 
is  P(G)  ^  for  all  4).  Analyzing  the  eigenvalues  of  G  for  each  error  mode  is  not  practical 
computationally.  Useful  approximations  are  adopted  here  to  simplify  the  problem.  We 
split  Z  =  Zj  +  ZR,  such  that  ^.(Zj)  are  gll  pure  imaginary,  and  X(ZR)  are  all  pure  real,  and 
analyze  the  individual  problems  G!  =  P^AtZj),  GR  =  PR(AtZR)  (where  now  Pj  and  PR  need 
not  necessarily  be  the  same  polynomial  function).  Considering  the  Fourier  matrix  Zj  (given 
in  Appendix  A),  we  estimate  : 

p(Gi)=  p(Pi(AtZi))  =  Pi(p(AtZj))  [3.24a] 

The  second  equality  in  equation  [3.24a]  is  strictly  true  only  if  AtZj  is  hermetian,  in  which 
case  [3.24a]  would  be  a  straightforward  application  of  the  spectral  mapping  theorem  (Varga 
(1962),  for  example).  In  the  present  context  of  local  linearized  analysis,  equation  [3.24a] 
is  submitted  as  an  approximation  which  provides  useful  information  as  will  be  seen. 

Matrix  ZR  is  pure  real  and  there  exists  S  such  that  Zr  =  S-1ZrS  is  real  and 
symmetric  (Abarbanel  and  Gottlieb  (1981)).  That  is,  the  spectral  mapping  theorem  applies 
directly  : 


P(Gr)  =  pCS^GrS)  =  p(PR(AtZR))  =  Pr(p(AiZr))  =  Pr(p(AiZr))  [3.24b] 

So  we  have  P(^r)  =  pR(p(AtZR))>  estimate  P(Gi)  *  Pj(p(AtZj))  It  remains 
then  only  to  ensure  that  the  locus  of  P(AtZi)  an(j  p(AtZR)  remajn  simultaneously  within  the 
scalar  stability  bound  for  the  chosen  scheme.  Accordingly,  if  we  evaluate  all  elements  of  Zj 


and  ZR  at  the  phase  for  which  they  are  a  maximum,  conservative  criteria  to  ensure  stability 
become : 
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Atp(Zi)  <  i  CFLop  and  Atp(ZR)<Qop 


[3.25] 


where  CFLop  and  £2op  are  operational  Courant  and  VonNeumann  numbers  chosen  such 
that  the  rectangle  formed  by,  0  <  Q.  <  Qop,  0  <  CFL  <  CFLop  lies  within  the  scalar 
stability  boundary  of  the  chosen  scheme  (see  Figure  3.2). 


So  it  is  assumed  in  equation  [3.25]  that  all  of  the  eigenvalues  of 

2  _  iS^Ai+iST1A2+iS^A3+i(-iDi) 


[3.26] 


are  imaginary,  and  that  all  of  the  eigenvalues  of 

ZR=2(Crl)S1+2(Q1-l)S2+2(C;-l)S3-S^S11T1-S4S;T2-ST1S;f3 


+4(Crl)2o4I+4(C11-l)2oT1I+4(Cel)2a;I+DR 


[3.27] 


are  real.  This  assumption  will  be  verified  below.  Here,  for  convenience,  D  has  been  split 

^  /S  0'S 

D  =  Di  +  Dr,  such  that  X(ZR)  are  all  pure  real  and  X(Zj)  are  all  pure  imaginary.  Di  and  Dr 
are  given  in  Appendix  A. 


Commencing  with  the  hyperbolic  analysis,  an  expression  for  Zj  (also  given  in 
Appendix  A)  can  be  obtained  from  equation  [3.26].  The  seven  eigenvalues  of  Zj  are : 
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X(Z)  =  iU, 


ju  ±^V 4o)2+4ccoki+c2(ki+k2+k3)±i-V 4co2-4co)k i  +c2{k  2 +k2+k2 ) , 


{u  *  kV(^+1|c,f,#+C2f2[ 


[3.28] 


where  U  s  US^+VS^+WS;  and  ki  s  ^xS^+ruS^+CxS^,  ki  ■  ^S^+riyS^+CyS^ 
k3  =  ^zSt+T|zSri+^zS)' 

**  1  These  eigenvalues  are  indeed  pure  imaginary  (as  per  our  choice  of 

Di),  since  the  groups  of  terms  inside  the  square  root  signs  are  never  negative  (independent 
of  the  sign  of  to).  The  coefficients  of  the  first  derivative  Jacobians  are  maximum  at  <J>  = 
jt/2.  A  conservative  estimate  for  the  maximum  value  of  p(Zj)  is  therefore  obtained  by 
evaluating  the  trigonometric  terms  at  this  phase,  so  the  spectral  radius  of  Zj  , 


p(Zj)  <  i  MAX 


{  U  +W 4co244c(oki-K:2(ki+k2+k3)+^V 4co2-4ccoki+c2(ki+k2+k3) 


u±fV(^+1Xc,f^-C2f2)J) 


[3.29] 


where  now  U  s  HWHW,  ki  =  fcNnJ+fcJ,  k2  =  kJ+ftsHcJ.  k3  =  feJ+hzHcJ. 


Moving  to  the  parabolic  stability  analysis,  we  have  : 


p(ZR)  =  p(  2 (C4-1)S,  +  2(Cn-I)S2  +  2(C?-1)S3 


-  S^S^Tj  -  S^T2-  Sr,S;T3 
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+  4(Crl)2a4I  +  4(CT)-1)2Ot1I  +  4(Crl)2^I  +  Dr  } 


[3.30] 


The  non-symmetric  element  of  Si  does  not  affect  the  eigenvalues  of  Si,  nor  of  ZR, 

«"S 

therefore  this  term  can  be  dropped  in  the  following  analysis,  so  Si  is  treated  as  symmetric 
and  real  hereafter.  We  can  invoke  the  matrix  inequality  : 


p(ZR)  <  ||  2(Clrl)Si  +  2(Ct1-1)S2  +  2(C^-1)§3 


-  S^S^Ti  -  S4S;T2-  Sr(S^T3 
+  4(q-l)2o4I  +  4(Cn-l)2o71I  +  4(Crl)2ocI  +  Dr  „ 


[3.31] 


where  1 1  ‘  II  represents  any  consistent  matrix  norm  (Stewart  (1973),  for  example).  The 
triangle  inequality  provides : 


p(ZR)  <  H2(Crl)Sill  +  II2(Cq- 1  )S2II  +II2(Cc-1)S3H 


+  ll-S^Tjll  +  ll-S^S^T2II  +  II-St,ScT3II 
+  l!4(q-l)2o4III  +  U4(Cn-l)2crT1III  +  H4(Crl)2a?III  +  IIDrII 


[3.32] 


The  coefficients  of  the  second  and  fourth  derivative  Jacobians  are  maximum  at 
<|>  =  it,  and  those  of  the  cross  derivative  terms  are  maximum  at  4>  =  tt/2.  A  conservative 
estimate  for  the  maximum  value  of  p(ZR)  is  therefore  obtained  by  evaluating  all  the 
trigonometric  terms  at  the  phase  for  which  they  are  the  largest : 
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p(ZR)  £  4llSill  +  4IIS2II  +  4IIS3II  +  IITill  +  IIT2II+IIT3II 


+  16lla^II!  +  1611^111  +  I6II05III  +  IIDrII 


Choosing  a  matrix  2  norm 


p(ZR)  <  4 [p(Si )  +  p(S2)  +  p(S3)]  +  [p(Tj)+  p(T2)+  p(T3)] 


[3.33] 


16[p(er$I)+  p(a11I)+p(a;I)]  +  p(DR) 


[3.34] 


where  we  have  used  the  matrix  identity  i  M I2  -  [p(wTw)]  ^  ^ 

/v  /\ 

Si,  Ti,  oj  and  Dr  ^  symmetric  and  real. 

✓V 

The  eigenvalues  of  Ti  are  complicated,  but  if  it  is  assumed  that  the  grid  is  close  to 
orthogonal  in  regions  where  the  viscous  stability  bound  becomes  important  (see  Martinelli 
(1987)),  the  eigenvalues  of  Tj  are  found  to  be: 

*i(Ti)  =  X2(Ti)  =  X3(Ti)  =  Um  =  XS(T  1)  =  0, 


U(Ji)  =  -XtCTO  -SV(V^VO(Vn*VT|) 


The  spectral  radius  is  therefore, 


[p(Tj)+  p(T2)+  p(T3)]  ”S£jV(v$.V£)(VtvVti) 


[3.35] 


'(V^v^)(vc.vo+V(VTiVn)(VC-vc)] 


[3.36] 
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A  simple  form  for  the  eigenvalues  of  S  i  can  be  obtained  without  invocation  of  a 

✓s 

local  grid  orthogonality  assumption,  but  for  consistency  with  the  development  for  Ti,  such 
terms  are  also  dropped  in  Si,  so  : 

Xi(Si)  =  0,  X2(Si)«4j2(V$-V$)f 

3p 


h(S,)  =  X4(Si) = ^<V5-  VU  Xj(S,)  -2^.+ vy , 
h(SD = ijm + £](V5-  vy,  hts,)  =  ijn, + vy 


For  y  =  1.4,  Prj  =  0.72,  Prt  =  0.9,  Prk  =  1.0,  Pre=1.3,  the  spectral  radius  is  therefore. 


[p(S, ,  +  p(S2)  +  P<S3>]-£(£+£M  +  ^  +  ^ 


[3.38] 


On  inspection : 

P(OtjI)+P(o^I)]  _  +  o^)  (recall  that  the  artificial  dissipation  scale 

factors  are  negative  and  real).  Lastly, 


x(dr)  =0, 0, 0, 0, 0, 2E  -  %  '  2C/iE  -  QU-rn 


Pk  p/: 


r  k 


pr 


[3.39] 


So 


p(dr)=maxJ 


2P  2Pi| 

Pk  pi2' 


2C?f?£  +  _^2e-yyj) 


[3.40] 
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The  spectral  radius  of  ZR : 


In  practice,  a  stable  local  timestep  could  be  obtained  from  : 


At  =  MIN  [  Ate ,  Atv  ]  s  MIN  [ 


v^r  i~>op  **op 

"pCZ^’p™ 


[3.42] 


where  P^l)  and  P(^*0  are  given  in  equations  [3.29]  and  [3.41].  It  is  noted  that  in  the 
absence  of  rotation  source  terms,  turbulence  model  source  terms,  and  artificial  dissipation 
terms,  expressions  [3.29]  and  [3.41]  represent  a  three-dimensional  extension  to  a  form 
presented  by  Martinelli  (1987). 


In  order  to  ascertain  the  relative  importance  of  the  various  terms  which  appear  in 
equations  [3.29]  and  [3.41],  an  order  of  magnitude  analysis  is  undertaken  here.  It  is  first 
observed,  upon  inspection  of  equations  [3.29]  and  [3.41],  that  the  direct  influence  of 
rotation  and  turbulence  source  terms  can  only  reduce  the  maximum  allowable  hyperbolic 
and  parabolic  timesteps,  for  plausible  values  of  co,  k ,  e  (-«>  <  co  <  «,  k  >  0,  e  >  0).  For 


the  following  arguments,  we  choose  constant  reference  length  scale,  /„,  =  blade  height  and 
velocity  scale,  =  blade  tip  speed 
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Upon  nondimensionalization,  the  terms  in  the  square  root  sign  on  the  second  line 
of  equation  [3.29]  are  of  order 


4Ro2  ±  ^Ni  +  jL(N?  +  Nl  +  N§) 


[3.43] 


_  coL  w  r  XT  L 

Ro  -  -y~,  M  -  N,  =  — — 

where  00  “  Axi  are  rotation,  Mach  and  grid  clustering  parameters.  By 

prescription  of  the  reference  scales,  V„  and  O(104)  <  R,,  <  0(10°).  For  the 
compressible  flows  considered  herein,  0(1(H)  <M  <  0(10°).  Notice  that  here,  M  does 
not  represent  a  typical  local  Mach  number,  rather  it  is  formed  from  a  reference  velocity  and 
local  speed  of  sound  For  realistic  grids,  O(10l)  <  N;<  O(K)5).  The  rotation  terms  will  be 
most  significant  when  R,,  is  large,  M  is  large  and  Nj  is  small,  that  is  in  the  core  flow 
region.  So  the  third  term  in  equation  [3.43]  is  no  less  than  one  order  of  magnitude  greater 
than  the  other  two  terms,  which  suggests  that  the  rotation  terms  can  be  neglected  in  the 
specification  of  a  local  timestep. 


To  substantiate  this  conclusion,  three  fully  turbulent  coarse  grid  test  rotor  flow 
cases  were  run.  The  nondimensional  stability  parameters  as  defined  above  are  Ro  =  0.5,  M 
s  0.2,  Nimin  =  10,  for  this  case.  In  light  of  the  above  arguments,  there  should  be  very 
litde  influence  of  the  rotation  terms  in  equation  [3.29]  on  the  stability  of  the  scheme.  In 
Figure  3.7,  convergence  rates  are  compared  for  the  cases  where  to  was  set  to  the  machine 
rotation  rate  of  1 13.2  s1  in  equation  [3.29],  and  where  co  was  set  to  zero  in  this  equation. 
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Figure  3.7  Comparison  of  convergence  rates  for  rotating  (solid)  and  stationary 
(dashed)  rotor  passage  flow  with  and  without  inclusion  of  rotation  terms  in  local 

dmestep  evaluations. 


The  plots  are  indistinguishable,  indicating  that  local  timesteps  are  indeed  negligibly 
influenced  by  system  rotation  for  this  case.  Accordingly,  the  first  term  in  equation  [3.29] 
may  be  replaced  by  the  expression : 
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'V^V^+Vti-Vti+VCVC+2[V^Vti+V^V;+VtiV;] 


[3.44] 


The  term  U  should  be  retained  in  equation  [3.44]  as  it  can  be  easily  shown  that  this 
term  may  be  of  the  same  order  of  magnitude  as  the  other  term  in  equation  [3.44]  if 
Miocal  =  q/c  =  1  Also  appearing  in  Figure  3.7  is  the  convergence  history  for  the  same  case 
run  with  the  machine  rotation  rate  set  to  zero  in  the  flow  solution.  The  convergence  rates 
compare  closely  to  approximately  a  five  decade  drop  in  the  mean  density  residual,  after 
which  the  rotating  case  converges  more  quickly.  This  result  is  included  to  stress  that  the 
physics  of  the  flowfield  can,  of  course,  influence  the  stability  of  the  scheme,  and  that  the 
above  conclusions  apply  only  to  the  determination  of  a  stable  local  timestep.  It  is  noted, 
that  these  conclusions  should  break  down  for  very  high  tip  speeds  and  very  coarse  grids,  in 
which  case,  for  stability,  the  rotation  terms  should  be  retained  in  equation  [3.29],  and  the 
rotation  source  terms  should  be  evaluated  at  every  stage. 


The  term  P/pe,  the  ratio  of  turbulent  kinetic  energy  production  to  dissipation, 
appears  in  equations  [3.29]  and  [3.41].  Assuming  that  the  flow  solver  converges  to  a 
physically  reasonable  solution,  it  can  be  argued  that  at  convergence  this  term  is  of  order  1 
or  less.  Specifically,  the  local  equilibrium  assumption,  P/pe  =  1,  is  valid  for  the  nearly 
homogeneous  core  flow  turbulence,  and  for  thin  shear  layers,  whereas,  near  the  boundary 
layer  edge  and  in  separated  shear  layers  and  wakes,  P/pe  <  1  (see  Hinze  (1975),  for 
example).  Additionally,  in  the  immediate  vicinity  of  the  wall  P/pCg  goes  to  zero,  where 
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subscript  "a"  denotes  the  true  dissipation  rate.  In  Chien’s  model,  the  dissipation  variable, 
e,  is  the  isotropic  component  of  the  dissipation  rate.  Though  the  isotropic  component  of 
the  dissipation  rate  goes  to  zero  at  the  wall,  unlike  the  true  dissipation  which  assumes  a 
finite  value  there,  P/pe  still  goes  to  zero.  To  see  this,  designate  subscripts  a  and  i  to 
correspond  to  actual  and  isotropic  components  of  the  dissipation  rate.  Then  : 

P  P 

P£i  P£a  +  D  [3.45] 

D  = 

where  /2  .  Both  numerator  and  denominator  in  equation  [3.45]  go  to  zero  as  the 

wall  is  approached.  The  denominator  going  to  zero  is  simply  a  statement  that  dissipation 
and  molecular  diffusion  of  turbulence  energy  are  in  balance  at  the  wall.  Through  the  entire 
near  wall  region  (say  y*  <  O(K)2)),  convection  of  turbulence  energy  is  negligible  by 
comparison  to  dissipation,  diffusion  and  production  (Hinze  (1975)).  Therefore,  in  order 
that  the  turbulence  energy  budget  be  in  balance  in  the  immediate  vicinity  of  the  wall, 

P 

Pe»  +  ® ,  must  uniformly  approach  zero  as  /  ~>  0.  So  with  Cifi  =  1-35,  1.4  <  Cih  ^1-8 
(for  Chien's  model),  it  can  be  assumed  that  at  convergence,  the  terms 


[3.46] 


in  equation  [3.29].  Turning  to  the  factor  E/k,  the  magnitude  of  the  hyperbolic  terms  in 
equation  [3.29]  are  now  compared.  In  the  core  flow  and  outer  regions  of  the  boundary 
layer 


L 

V„ 


cVv^v^+vti-vti+vc-v;+2[V^vti+v^vc+vti-v;]  =  ojJ^  n»)  >  101 


[3.47] 
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where  the  same  constant  reference  length  and  velocity  scales  are  chosen  as  before.  If  local 
turbulence  intensity  and  turbulence  length  scales  are  defined : 


T.e.i.s£ 

e 


[3.48] 


then  V  00 '  ,  where  C  is  the  ratio  of  integral  turbulence  length  scale  to  blade 

height, L  =  a~.  Reasonable  bounds  are  O(10  3)  <  T  <  O(10-2),  O(l0'3)  <  C  <  O(10-2). 

[iik  <io° 

So  v  Voo/K'  ,  that  is,  comparing  equation  [3.47],  awav  from  the  immediate  vicinity  of 
solid  boundaries,  the  turbulence  source  terms  should  be  no  more  than  one  order  of 
magnitude  less  than  the  convective  acceleration  terms  in  equation  [3.29],  at  convergence. 


In  the  immediate  vicinity  of  the  blade  and  endwalls,  an  appropriate  near  wall  length 
scale,  v/u*,  and  velocity  scale,  u*,  are  chosen.  The  ratio  of  convection  to  turbulence 
source  terms  is  of  order 

c/Ay  _  (c/u*)  ([v/u*]/Ay) 

e/k  (e+/k+)  [3.49] 

where  y  is  chosen  as  the  blade  normal  coordinate  here  for  familiarity, u  =  is  the 
e+  =  k+  =  j£_ 

friction  velocity,  and  u*4,  u*2.  Again,  in  Chien’s  model,  £  is  the  isotropic 

component  of  dissipation  rate.  For  turbulence  models  which  solve  for  the  actual 
dissipation  rate  (Lam-Bremhorst  (1981),  Myong  and  Kasagi  (1990)  for  instance), 

£+/k+  — >  00  as  y  — >  0,  since  e  is  finite  at  a  solid  boundary.  This  may  be  why  these  models  show 
unfavorable  stability  properties  in  time-marching  Navier-Stokes  computations,  as  noted  by 
the  present  author  and  others  (Choi  and  Knight  (1989),  Kirtley  (1989)).  In  the  opinion  of 
the  author,  this  near  wall  behavior  is  more  likely  the  cause  of  numerical  problems  than  the 
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Neumann  or  second  derivative  boundary  conditions  on  the  e  equation,  to  which  these 
difficulties  are  often  attributed,  in  models  which  solve  for  e a  (Speziale  et  al.  (1990),  for 
example). 

For  Chien’s  model,  we  estimate  from  published  near  wall  distributions  of  k+  and 
e+  for  simple  shear  layers  (Patel  et  al.  (1985))  the  following : 


Also 


0(1 02)  <  E+/k+  <  0(1 01) 


OdO"1)  <  ([v/u*]/Ay)  *  ([y/Ay]/y+)  <  0(10°) 


[3.50] 


[3.51] 


depending  on  near  wall  grid  resolution  (at  the  wall  ([y/Ay]/y+)  =  l/y+  =  0(10°) ),  and 


(c/u*)  =  (c/Vo.) 


vyr 


B(#|- 


OdO1) 


[3.52] 


So  near  the  wall,  we  also  expect  the  turbulence  source  terms  to  be  no  more  than 
one  order  of  magnitude  less  than  the  convective  acceleration  terms  in  equation  [3.29),  at 
convergence. 


Turning  to  the  other  turbulence  source  terms,  we  can  invoke  exactly  the  same 

2E(£j  2q  f  [gl 

arguments  as  above  to  neglect  the  terms  Pe  ^ '  and  2  'k  >  which  appear  in  equation 
[3.41],  since  they  will  be  negligible  by  comparison  to  the  convective  acceleration  terms  in 


equation  [3.29]  everywhere  in  the  flowfield.  However,  the  near-wall  modelling  terms, 

2[Xic.y.„ 

l 2  and  / 2 


can  be  very  large  near  solid  boundaries  and  should  be  considered.  Again 
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choosing  y  as  the  blade  normal  coordinate,  near  solid  surfaces  the  ratio  of  physical 

Hi  +  Ht  /  Hi 

dissipation  to  wall  correction  terms  in  equation  [3.41]  is  of  order  Ay2  y2  'phis 
expression  will  be  very  large  except  in  the  immediate  vicinity  of  the  wall  ((^  — >  0,  y  --> 

Ay),  where  it  approaches  unity.  Therefore,  very  near  a  solid  wall,  where  turbulent 
diffusivity  damping  effects  dominate,  the  wall  modelling  functions  give  rise  to  parabolic 
stability  restrictions  approximately  equal  in  magnitude  to  physical  diffusion  terms  (This  is 
not  a  surprising  result  since  D  represents  molecular  diffusion  of  k  at  the  wall).  However,  it 
is  shown  in  the  next  section  that  convection  dominates  the  overall  stability  of  the  scheme  in 
such  regions,  so  it  is  expected  that  the  wall  modelling  terms  in  equation  [3.41]  may  be 
neglected  as  well. 


In  summary,  the  foregoing  arguments  suggest  all  turbulence  source  terms  have 
negligible  influence  on  the  numerical  stability  of  the  scheme,  at  convergence.  To 
investigate  this,  the  relative  magnitudes  of  the  terms  above  were  examined  for  the  fine  grid 
rotor  and  cascade  calculations  in  both  early  and  late  stages  of  iteration.  Designating 


A) 

MpMCif,£+C2f2) 

B) 

2 E£ 
pe  k 

2pi 

C) 

2Czf2k 

D) 

P  l2 

E)  U  +  cVv^- V^+Vt)-Vti+VC-VC+2[V^Vti+V^V;+Vti- VC] 


5(a^)n.ntv,v,+vt.vt) 

+  ^■[V(V5V5)(VT)Vn)W(V5V5)(v;v;)W(VnVTi)(vcvo] 
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For  the  fine  grid  rotor  flow  computation,  at  convergence,  the  magnitudes  of  terms 
A,  B  and  C  were  less  than  0.1  times  the  magnitude  of  term  E,  and  the  magnitude  of  term  D 
was  less  than  0.1  times  the  magnitude  of  term  F,  at  every  point  in  the  computational 
domain,  thereby  substantiating  the  arguments  of  the  previous  section  for  this  case.  Upon 
cold  start  "inviscid"  initialization  the  same  conditions  were  observed  except  for  term  B. 
Specifically,  at  iteration  1,  the  magnitude  of  term  B  was  greater  than  one  tenth  the  value  of 
term  E  at  approximately  .15  of  the  grid  points  in  the  flowfield.  This  number  decreased  to 
.002  in  250  iterations.  For  the  cascade  computation,  the  ratio  of  these  terms  were  all  less 
than  0.1,  both  at  iteration  1  and  at  convergence. 

Some  discussion  is  warranted  here.  The  foregoing  arguments  suggest  that  the 
stiffness  associated  with  large  values  of  near  wall  turbulence  transport  source  terms  is 
negligible  for  Chien's  model,  except  in  early  stages  of  iteration.  The  stiffness  (local 
timestep  restriction)  associated  with  computing  flows  on  very  highly  stretched  grids  using 
time  marching  algorithms  supersedes  difficulties  brought  about  by  these  source  terms.  For 
this  reason,  the  present  author  has  been  able  to  compute  successfully  a  variety  of  high 
Reynolds  number  flowfields  using  this  model  (see  following  chapters),  without  taking 
account  of  the  k-e  source  terms  in  the  determination  of  a  local  timestep  and  while  evaluating 
source  terms  prior  to  the  first  stage  only.  Additionally,  the  convergence  rates  have  been 
found  to  be  nearly  identical  to  calculations  performed  using  an  algebraic  eddy  viscosity 
model.  This  is  illustrated  in  Figure  3.8,  which  has  been  adapted  from  Kunz  and 
Lakshminarayana  (1990).  In  this  figure,  a  flat  plate  plate  boundary  layer  computation  was 
performed  on  a  highly  stretched  grid  (refer  to  validation  studies.  Chapter  4)  using  both  the 
low  Reynolds  number  k-e  model  and  an  algebraic  eddy  viscosity  model  due  to  Baldwin  and 
Lomax  (1978).  As  seen  in  the  figure,  the  convergence  rates  are  very  similar. 


Iteration# 


Figure  3.8  Convergence  histories  for  turbulent  flat  plate  flow  computations. 
Baldwin  and  Lomax  (1978)  model  (solid  line)  and  present  k-e  model  (dashed 
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Since  the  source  terms  do  not  significantly  influence  the  stability  of  the  scheme,  the 
use  of  locally  implicit  treatment  of  the  source  terms  should  not  be  necessary.  This 
reconciles  with  recent  results  obtained  by  Martelli  and  Michelassi  (1990)  who  used  a 
factored  implicit  time  marching  scheme.  They  found  that  pointwise  implicit  treatment  for 
the  source  terms  in  the  near-wall  q-co  model  they  used  provided  no  difference  in 
convergence  rate  over  a  purely  explicit  treatment. 

To  investigate  this  further,  the  fine  grid  cascade  computation  was  run  three  times. 
First  the  standard  approach  was  used,  where  all  source  and  dissipative  terms  are  evaluated 
prior  to  the  first  stage.  Next,  the  code  was  run  evaluating  all  source  and  dissipative  terms 
at  every  stage.  Lastly,  the  code  was  run  using  pointwise  implicit  treatment  of  the 
turbulence  source  terms.  Specifically,  for  the  k  and  e  equations,  each  stage  in  the  standard 
RK4  procedure  (equation  [3.4])  is  replaced  by  : 

Qk+1  =  Q°  +  ock+1At  [i  -  ctk+iAtDk] 1  R(Qk),  k  =  0,  N-l  [3.53] 

where  D  is  given  in  Appendix  A  (note :  0)  =  0  for  this  case).  For  each  of  the  three 
computations,  CFL^  was  chosen  as  2.8,  and  £2op  was  adjusted  to  the  maximum  value  for 
which  a  stable  solution  could  be  obtained.  Specifically,  =  0.95,  -  0.80, 

£2op3  =  0.30.  The  convergence  histories  for  these  calculations  appear  in  Figure  3.9a.  In 
Figure  3.9b,  the  predicted  boundary  layer  velocity  and  turbulent  kinetic  energy  profiles  at 
s/c  =  0.2  on  the  suction  surface  are  compared.  The  convergence  histories  are  quite  similar 
for  all  three  approaches,  and  the  converged  solutions  are  indistinguishable.  These  findings 
further  substantiate  the  claims  made  above. 

Recently,  Gerolymos  (1990)  obtained  two  dimensional  supersonic  duct  flow 
results  on  very  fine  grids  using  an  explicit  multigrid  scheme  and  a  low  Reynolds  number 


#  iterations 


V/Ved*  V  Vj* 


Figure  3.9  Comparison  of  a)  convergence  histories,  b)  converged  boundary  layer 
velocity  and  local  turbulence  intensity  profiles  at  a  location  on  the  suction  surface  far 
fine  grid  cascade  computations.  Results  for  three  different  source  and  dissipation  term 

treatments  are  compared. 
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model  due  to  Launder  and  Sharma  (1974).  This  model  is  similar  to  Chien’s  in  that  the 
isotropic  component  of  dissipation  is  a  transport  variable.  Gerolymos  notes  that  "for  both 
the  driver  scheme  and  the  multiple-grid  procedure,  [pointwise  implicit  treatment]  of  the 
source  terms  is  essential  to  the  stability  of  the  method."  Since  the  source  term  linearization 
he  chose  is  similar  to  that  analyzed  and  used  above,  his  statement  for  tne  driver  scheme  is 
in  contradiction  to  the  present  findings.  However,  when  a  multigrid  procedure  is 
incorporated,  the  thesis  conclusions  reached  above  become  less  valid  if  the  complete 
residuals,  including  source  terms,  are  evaluated  on  successively  coarser  grids. 

Specifically,  grid  metric  stability  constraints  associated  with  convection  and  diffusion 
operators  are  relaxed  when  going  to  coarser  grids,  whereas  timestep  constraints  associated 
with  the  source  terms  are  not.  This  is  pointed  out  by  Mavriplis  and  Martinelli  (1991),  who, 
like  Gerolymos,  chose  a  multigrid  procedure  where  the  turbulence  source  term  residuals  are 
retained  on  coarse  grids.  They  justify  using  implicit  source  treatment  so  that  these  source 
terms  will  not  limit  the  larger  timesteps  allowed  by  relaxing  the  convection  and  diffusion 
stability  constraints  on  coarser  grids.  This  exact  same  argument  may  be  extended  to 
consideration  of  the  rotation  source  terms  as  well. 

In  practice  the  author  has  encountered  several  cases,  including  flat  plate  and 
supersonic  cascade  flow  computations  where  large  relative  values  of  production,  P/pe, 
cause  the  turbulence  model  to  become  rapidly  unstable  in  early  stages  of  iteration,  when 
standard  inviscid  initialization  is  used,  These  cases  have  in  common  that  freestream 
turbulence  levels  are  low  and  Reynolds  numbers  are  very  high.  It  has  been  found  that 
incoxporation  of  the  turbulence  source  terms  in  equations  [3.29]  and  [3.41]  is  ineffective  in 
stabilizing  these  cases.  However,  several  procedures  are  effective.  Specifically,  k  and  e 
are  required  to  remain  positive : 


pe  >  K£  p.oE.0 


[3.54a] 
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pk  >  Kk  p  Jc.  [3.54b] 

where  Ke  =  0.01  -  0.0001  and  Kk  =  0.0001.  These  two  procedures  were  recommended  to 
the  author  by  Connell  (1989),  and  are  used  by  Gerolymos  (1990).  Also  the  inflow 
turbulence  intensity  can  be  increased  to  an  unrealistically  high  value  for  the  first  20-50 
iterations : 


T  =  5 


actual 


[3.54c] 


This  has  the  effect  of  significantly  reducing  P/pe  in  early  stages  of  iteration,  since  upon 


initialization 


-E-clKil-o 

pe  '  e/ 


'll' 

A 


So  in  lieu  of  the  above  analyses  and  verification  studies,  rotation  and  turbulence 
source  terms  are  not  included  in  the  specification  of  a  local  timestep,  in  any  of  the 
computations  to  follow. 


3.3.3.3Phvsical  and  Artificial  Dissipation 

Turning  next  to  the  physical  diffusion  terms  in  equation  [3.41],  both  of  order 


,  where 


pVJoe 

Rel*1-—  ,Rei 
Hi 


pvj L 

lit 


[3.55] 


In  many  internal  and  external  flows  including  turbomachinery  flows,  we  have 

Rei »  C)(l05  -  106)  ^  (ioi)  <  <  OflO5),  except  in  regions  where  wall  damping  is 

effective  (where  Ret  -->  «>).  So  comparing  the  inviscid  terms  in  equation  [3.29],  which  are 
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JL  n- 

of  order  M  \  we  invoke  the  bounds  on  M  and  N;  provided  in  the  previous  sections,  and 

conclude  that  the  parabolic  timestep  may  be  smaller  than  the  hyperbolic  term,  Atv  <  A^,  in 

regions  of  very  high  grid  clustering  (Nj  =  O(K)5))  and  low  (  =  CK101)), 

_L_  n?  >  _L  Ni 

ie.  Ret  1  M  .In  Figure  3.10,  the  Aty  /Atc  =  1  contour  is  plotted  for  the  fully 


converged  production  and  coarse  grid  rotor  computations  at  an  annular  grid  slice  near 
midspan.  In  both  cases,  the  parabolic  stability  considerations  described  above  are  seen  to 
manifest  themselves  near  the  leading  and  trailing  edges.  This  provides  numerical 
verification  of  the  necessity  to  include  the  effective  diffusivity  terms  in  equation  [3.42]  in 
calculations  where  large  diffusivities  occur  in  regions  of  large  grid  stretching.  It  has  been 
the  experience  of  the  author  that  this  is  especially  important  during  the  course  of  iteration  as 
illustrated  in  Figure  3.11.  There  it  is  shown  that  the  parabolic  stability  bound  dominates 
for  more  than  half  of  the  grid  points  in  early  stages  of  iteration  for  the  coarse  grid  test  flow 
configuration,  and  for  approximately  7%  of  the  grid  points  at  convergence.  The  cascade 
computation  yielded  similar  results.  Specifically,  at  convergence  At*  <  At,,  for 


approximately  4  %  of  the  grid  points,  all  in  the  immediate  vicinity  of  the  blade  leading  and 
trailing  edges.  It  is  noted  that  algebraic  eddy  viscosity  models  do  not  typically  give  rise  to 
the  parabolic  stability  mechanism  presented  here,  when  applied  to  turbomachinery  blade 
row  computations.  This  is  because  eddy  viscosity  is  assumed  to  be  zero  upstream  of  the 
blade  row  (as  for  say  Baldwin  and  Lomax  (1978)  point  transition  model),  so  that, 
depending  on  wake  treatment,  there  should  be  no  regions  of  the  flowfield  where  iarge  eddy 
viscosity  is  coincident  with  high  levels  of  grid  clustering. 


Artificial  dissipation  is  added  to  the  discrete  equations  by  design  to  influence  the 
stability  of  the  scheme.  If  physical  dissipation  gives  rise  to  a  locally  dominant  parabolic 
stability  constraint,  any  artificial  dissipation  will  destabilize  the  calculation.  Therefore,  the 
artificial  dissipation  terms  in  equation  [3.41]  must  be  retained  if  maximum  operational 


and 


b)co^Sridmtmflow ca)  ^  Wh“'  *.  fora,  ft* 

Widens  gnd  siice 


Figure  3. 1 1  Percentage  of  grid  points  for  which  At*  <  Atc  as  a  function  of 
iteration  number  for  coarse  grid  rotor  flow  computation. 
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VonNeumann  numbers  are  to  be  used  independent  of  levels  of  added  dissipation.  This  is 
illustrated  in  Figure  3. 12.  The  solid  line  in  this  figure  is  the  convergence  history  for  the 
coarse  grid  test  rotor  flow  computation  run  with  Qop  -  1.75,  CFLop  =  2.80.  Standard 
values  of  fourth  difference  dissipation  (refer  to  section  3.2),  were  used  for  this  calculation, 
the  local  timesteps  being  computed  from  equations  [3.29]  and  [3.41].  Due  to  the  parabolic 
stability  constraints  arising  from  physical  dissipation,  Qop  had  to  be  lowered  substantially 
when  artificial  dissipation  was  not  incorporated  in  equation  [3.41].  The  three  dashed 
curves  in  Figure  3.12  illustrate  this.  A  choice  of  Qop  =  .6  yielded  a  divergent  solution;  £lop 
=  .5  showed  some  unstable  behavior  in  early  stages  of  iteration  (compare  Figure  3.1 1),  but 
recovered  and  converged.  £20p  =  .4  provided  a  stable  iterative  process  with  a  slightly 
diminished  convergence  rate.  This  demonstrates  that  it  is  useful  to  incorporate  artificial 
dissipation  in  the  parabolic  timestep  prescription  in  problems  of  this  type,  to  avoid  ad  hoc 
specification  of  an  optimum  VonNeumann  number  which  varies  with  amount  of  artificial 
dissipation  used. 

3., 3.4...S, lability  Boundary 

The  investigations  reported  in  the  previous  section,  suggest  that  the  rotation  and 
turbulence  source  terms  can  be  neglected  in  equations  [3.29]  and  [3.41],  but  all  other  terms 
including  artificial  dissipation  should  be  included,  and  the  procedures  summarized  in 
equations  [3.54a-c]  should  be  adopted.  This  gives  rise  to  a  robust  and  convergent  code  for 
computation  of  internal  flows,  including  turbomachinery  flows,  using  the  low  Reynolds 
number  k-e  model  in  an  explicit  solution  procedure.  In  Figure  3.13,  the  location  of  CFL^ 
and  Qgp  corresponding  to  maximum  attainable  convergence  rates,  for  the  fine  grid  cascade 
and  rotor  flow  computations,  are  plotted  with  the  stability  boundary  for  the  scheme.  The 
locations  of  these  operating  points  are  very  close  to  the  scalar  bound  for  the  scheme 
(compare  Figure  3.2).  The  corresponding  convergence  histories  for  these  two 
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#  iterations 


Figure  3.12  Comparison  of  convergence  rates  for  several  operating  Von 
Neumann  numbers  with  (solid)  and  without  (dashed)  inclusion  of  artificial 
dissipation  terms  in  timestep  computations. 
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Figure  3.13  Location  of  numerical  operating  point,  (CFLop,  H^),  for  full  scale  cascade  and  rotor  flow 
computations,  plotted  in  the  complex  Z  plane  with  the  stability  boundary  for  the  RK4  scheme  used. 
Corresponding  convergence  histories  for  these  two  calculations  are  also  included. 
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"production"  runs  also  appear  there.  This  figure  provides  additional  support  to  the  validity 
of  the  foregoing  analyses  in  application  to  full  scale  engineering  computations. 

3.4  Numerical  Implementation  and  Stability  of  Algebraic  Reynolds  Stress  Model 

In  Chapter  2,  a  hybrid  algebraic  Reynolds  stress  model  was  introduced.  Equations 
[2.14]  are  implemented  directly  in  the  flow  solver  as  follows.  The  nonlinear  coupled 
algebraic  system  defined  by  equation  [2.14]  can  be  implemented  by  lagging  the  turbulence 
energy  production  term,  P,  by  one  iteration  or  by  using  the  production  of  turbulence  kinetic 
energy  obtained  from  the  2-equation  model.  The  latter  approach  is  adopted  here.  If  the 
machine  axis  is  chosen  to  be  coincident  with  the  x-axis  (no  loss  of  generality),  the  cartesian 
components  of  the  rotation  production  tensor  (equation  [2.14])  are  simplified.  The 
resulting  6x6  algebraic  linear  system  of  equations  become  A  R  =  B  where  : 
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A  = 


1+L[k2(|ux)] 

Lfe2(-|vy)] 

L[k2(-|wz)] 

L[k2(-|ux)] 

1+Uk2(|vy)] 

Uk2(-|wz)] 

L[k2(-|ux)] 

L[k2(-|vy)] 

1+L[k2(|wz)] 

L[k2(vx)] 

L[k2(uy)] 

0 

L[k2(wx)] 

0 

L[k2(uz)] 

0 

L[ki(2Q)f 

k2(wy)] 

L[ki(-2n>f 

k2(vz)] 

[3.56] 

U^uy^vx)]  L[k2(|ur|wx)]  L^2(-|vz-2-wy)] 
Uk2(|vxf  y)l  LIK^WJ] 

'**&$*  ^<fwx^)] 

1  +L[k2(v  y+Ux)]  Lk2(v^  Uk^Uz)] 

1+L[k2(Wz+Ux)]  Lfk2(uy)l 

Uk2(vx)]  1+L{k2(wz+Vy)] 


,  with  Cj  =  1.5,  C2  =  0.6  as 


R  =  [-pu'iu'i5  -pu'2u'2>  -pu'3u3>  -pu'iu2>  -puiu3>  -pu'2u3]T 


B  =  ['3pk,  '3pk,  "3pk,  0,  0,  0]T 


.  ’  ki=^fe,  k2=l-C2,  k2=Ci-l 
and  p  +  Pek3  2 


empirically  obtained  (Launder  et  al.  (1975))  FRS  modelling  constants. 


L[k2(wx)] 


System  [3.56]  is  solved  at  every  grid  point  directly  by  Gaussian  elimination 
[subroutine  provided  by  Rao  (1982)]. 

The  direct  inversion  of  a  linear  6x6  system  of  equations  at  every  grid  point,  at 
every  iteration  would  be  prohibitive  computationally.  Numerical  experimentation  has 


shown  that  evaluating  the  Reynolds  stresses  every  20  iterations  (while  reevaluating  viscous 
fluxes  every  iteration  for  stability)  gives  rise  to  a  convergent  procedure. 
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Recalling  the  effective  stress  tensor  for  the  ARS  model : 


xy  =  m 


2  s..d«k' 

,'dxj  dxj 

'  3  6,J  3xkJ 

-  pkTy 


[2.15] 


By  comparison,  the  effective  stress  tensor  for  the  two-equation  model  used  can  be  written  : 


So  for  the  ARSM,  an  "anisotropic  eddy  viscosity"  can  be  identified  : 


M-tij  = 


[3.58] 


It  has  been  found  that  this  identification  is  useful  in  evaluating  the  numerical 
stability  of  the  turbulence  model.  Specifically,  if  the  individual  components  of  as 
defined  by  equation  [3.58],  are  all  of  the  same  order  of  magnitude  as  (it  provided  by  the  k-e 
model,  then  the  effective  diffusivity  arguments  provided  in  the  stability  analysis  of  Section 
3.3.3.3  remain  valid.  That  is,  the  eddy  viscosity  computed  in  the  k-e  solution  can  still  be 
used  to  provide  a  local  stable  timestep,  and  the  convergence  rate  of  the  numerical  scheme 
should  not  be  significantly  affected,  when  the  ARSM  is  used. 


To  illustrate  this,  a  flat  plate  turbulent  flow  computation  was  performed  using 
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1)  the  low  Reynolds  number  k-e  model,  2)  the  pure  high  Reynolds  number  ARSM,  3)  the 
hybrid  model,  with  y+match  =  200.  The  convergence  history  for  the  computation  using  the 
hybrid  model  showed  some  "peaking"  behavior  during  the  course  of  iteration.  This 
phenomena  was  traced  to  grid  points  whose  local  y+  values  were  very  close  to  y+matrh 
during  iteration.  If  this  local  value  crosses  y+match*  the  other  turbulence  model  is  suddenly 
used  locally,  and  this  gives  rise  to  a  very  large  local  residual  which  in  turn  strongly  affects 
the  RMS  residual.  To  alleviate  this  problem,  a  weighting  function  approach  was  devised, 
which  blends  the  two  models  near  y+malch.  Specifically  the  viscous  residuals  are  weighted 
as  follows : 

/N  /\  A 

Rviscous  =  W  ARSM^-ARSM  +  W^-eRk-e 


where 


Warsm  =  2 


Wk-£  =  1.  -  WArsm 


[3.59] 


These  weighting  functions  are  plotted  in  Figure  3.14,  for  a  clustering  parameter  (3  =  10.0, 
which  is  the  value  used  in  the  calculations  presented  using  the  hybrid  model  in  this  thesis. 
In  Figure  3.15,  the  convergence  history  for  the  blended  hybrid  model  is  seen  to  closely 
compare  to  those  of  the  two  component  models,  thereby  substantiating  the  stability 
argument  presented  above  for  this  case. 

In  the  first  section  of  the  next  chapter,  the  results  of  these  flat  plate  computations 


are  presented. 


5000  10000 

Iteration  # 

Figure  3.15  Convergence  histories  for  turbulent  flat  plate  flow  computations. 
Low  Reynolds  number  k-e  model  (alternating  long  dash-short  dash),  high 
Reynolds  number  ARSM  model  (solid),  hybrid  model,  blending  function, 
equation  [3.63]  used  (short  dash). 
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At  the  time  of  the  writing  of  this  thesis,  the  vectorized  code  executed  at  4.2  x  lO"5 
CPU  seconds  /  (gridpoint  *  iteration)  on  the  NAS  Cray  Y-MP  at  NASA  Ames  Research 
Center,  when  the  low  Reynolds  number  k-e  model  is  used.  As  discussed  above,  the 
convergence  penalty  for  implementing  a  two-equation  model  in  a  purely  explicit  procedure 
is  small.  The  two  additional  transport  equations,  which  require  only  partially  vectorizable 
source  term  and  eddy  viscosity  computations  give  rise  to  an  overall  factor  of  1.5  in 


execution  rate  over  laminar  flow  computations  on  a  Cray  2.  The  current  overhead 


associated  with  computations  using  hybrid  k-e/ARS  model  is  an  additional  40  %  for 


three-dimensional  calculations  on  the  Y-MP.  This  overhead  is  due  to  computation  of  both 
k-e  and  ARSM  viscous  flux  residuals  at  every  iteration,  and  re-evaluation  of  the  Reynolds 


stress  components  every  20  iterations. 
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CHAPTER  4 

PREDICTION  OF  TWO-DTMENS TONAL  INTERNAL  AND  CASCADE  FLOWS 

In  this  chapter,  results  and  interpretation  of  two-dimensional  flow  computations 
are  presented. 

Two  complex  cascade  flows  are  computed,  for  supersonic  and  low  subsonic 
freestream  conditions.  For  the  supersonic  cascade,  isentropic  blade  Mach  number,  shock- 
boundary  layer  structure  and  wake  loss  profiles  are  compared  with  experimentally 
measured  values.  Pressure  distribution  and  boundary  layer  profiles  of  velocity  and 
turbulent  kinetic  energy  are  compared  with  data  for  the  subsonic  cascade. 

4.1  Turbulent  Flat  Plate  Boundary  Laver  Flow  -  Turbulence  Model  Validation 

In  Chapters  2  and  3,  the  formulation  and  numerical  implementation  of  the 
turbulence  models  incorporated  in  this  thesis  were  presented.  Stability  analysis  and 
numerical  verification  results,  including  convergence  histories  for  simple  shear  layer 
computations,  using  a  compressible  extension  of  Chien's  (1982)  low  Reynolds  number  k-e 
model  and  a  hybrid  k-e/ARS  model,  were  presented.  It  remains  to  present  accuracy 
validation  for  these  computations.  These  results  are  provided  in  this  section. 

A  33  x  33  computational  grid  was  used  for  the  flat  plate  computations.  The  grid  and 
freestream  conditions  were  chosen  to  provide  a  Reynolds  number  based  on  boundary  layer 
thickness,  Rgg  =  80000,  near  the  exit  of  the  computational  domain.  This  corresponds  to 
the  flat  plate  measurements  of  Klebanoff  (1954),  against  which  the  computational  results 
provided  here  are  compared. 


100 


For  both  turbulence  models  used,  a  low  Reynolds  number  k-e  model  is  used  in  the 
very  near  wall  region.  Accordingly,  it  is  appropriate,  where  practically  feasible,  to  provide 
near  wall  grid  clustering  such  that  y+  =  1  at  the  first  grid  point  off  of  solid  boundaries. 

This  reconciles  with  the  observations  of  many  authors  (see  Awa  et  al.  (1990),  for 
example),  that  grid  independent  solutions  using  such  models  require  such  near  wall  grid 
spacing. 

To  ensure  that  the  first  grid  point  off  of  solid  boundaries,  will  be  near  a 
prespecified  value  of  y+,  a  standard  empirical  momentum  integral  relation  for  skin  friction 
on  a  turbulent  developing  flat  plate  flow  is  used  (see  Fox  and  McDonald  (1985),  for 
example) : 


Twall  ~ 


^£(0.0577) 

(Rex)1/5 


[4.1] 


AyWai,V^f 

Using  the  definition  of  y+,  ysPecified  v  ,  an  expression  for  Aywau  can  be 

obtained : 


Ay  wall  —  specific 


(ReL)1/5 


\l/2 


^(0.0577) 


[4.2] 


„  V_L 
Rei  - - 

where  v  ,  and  L  is  taken  as  the  plate  length  (or  chord  length,  as  it  were).  Because 

AyWall  *s  a  weak  function  of  Reynolds  number  (proportional  to  Re^1/10),  equation  [4.2] 
provides  an  appropriate  value  of  Aywau  along  the  length  of  the  plate. 
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The  grid  spacing  at  the  wall  for  the  flat  plate  computation  was  3.1  x  lO-6  times  the 
plate  length.  The  aspect  ratio  of  near  wall  grid  elements  in  this  extremely  highly  stretched 
grid  reached  30000  at  the  aft  end  of  the  plate.  Uniform  total  pressure  and  temperature 
profiles  .vere  specified,  and  the  back  pressure  was  set  such  that  the  free  stream  Mach 
number  was  approximately  0.5  to  accommodate  the  compressible  formulation  of  the  flow 
solver.  Three  calculations  were  performed,  using  the  low  Reynolds  number  form  of  the 
k-e  model,  the  high  Reynolds  number  form  of  the  ARSM  with  no  near  wall  treatment,  and 
the  hybrid  model,  with  y+ma[ch  =  200.  As  shown  in  Chapter  3,  the  artificial  dissipation 
levels  were  scaled  by  the  square  of  the  ratio  of  local  velocity  to  free  stream  velocity. 

In  Figure  4.1,  the  predicted  velocity  profiles  are  plotted  in  inner  variables  at 
Reg  =  80000,  and  compared  against  standard  law-of-the-wall  correlations.  The  k-e  and 
hybrid  models  match  the  correlations,  whereas  the  high  Reynolds  number  ARSM  shows 
poor  agreement  as  expected.  Klebanoff  (1954)  measured  the  turbulent  kinetic  energy 
profile  in  the  boundary  layer,  and  his  data  are  used  for  comparison  with  predictions  of  this 
parameter  in  Figure  4.2.  Again  agreement  is  excellent  for  the  k-e  and  hybrid  models. 

In  Figure  4.3,  the  distribution  of  turbulence  intensities,  predicted  using  the  hybrid 
model  are  shown,  and  compared  to  Klebanoff  s  measurements.  Comparison  is  good,  with 
’’e  notable  exception  that  the  hybrid  model  does  not  predict  the  distribution  of  energy 

•*  H  M  ** 

between "  PU2U2, "  PU3U3.  This  is  a  direct  consequence  of  not  explicitly  modelling  the 
pressure  strain  term.^'i-w  ,  presented  in  Chapter  2.  Despite  this  modelling  shortcoming,  the 
difference  between  the  streamwise  and  wall  normal  intensities,  is  reasonably  well 
predicted,  a  prerequisite  for  predicting  the  influence  of  rotation  on  the  turbulence  structure 
in  a  rotating  boundary  layer.  Specifically,  Reynolds  shear  stress  source  terms  in  the  full 
Reynolds  stress  and  algebraic  Reynolds  stress  equations  are  directly  proportional  to  both 
the  difference  between  these  intensities,  and  the  rotation  rate.  The  location  of  y^match ls 
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indicated  in  Figure  4.3.  Preliminary  efforts  to  validate  the  hybrid  model  for  an  isolated 
rotational  strain  case  were  undertaken.  Specifically,  a  rotating  fully  developed  channel 
flow,  tested  experimentally  by  Johnston  et  al.  (1972),  was  computed.  An  appropriate 
deflection  of  the  mean  flow  was  observed,  as  was  the  characteristic  redistribution  of 
Reynolds  intensities  near  the  leading  and  trailing  sides  of  the  channel.  Unfortunately  a 
converged  solution  for  this  flow  was  not  attainable  due  to  the  compressible  formulation  of 
the  flow  solver,  and  the  fully  developed  nature  of  the  measured  channel  flow.  It  is 
anticipated  that  appropriate  validation  against  a  developing  boundary  layer  subject  to 
spanwise  rotation  will  provide  more  precise  validation  in  regards  to  isolated  rotational 
strain. 


The  following  cascade  computations  were  performed  with  an  earlier  version  of  the 
code ,  where  the  k  and  e  equations  were  numerically  decoupled  (as  defined  in  Chapter  3) 
from  the  mean  flow  equations.  Additionally,  the  cases  were  run  in  laminar  mode  for 
several  hundred  iterations  before  the  k-e  model  was  turned  on  (see  convergence  histories 
below).  This  strategy  was  used  by  the  author  in  earlier  versions  of  the  code,  to  help  reduce 
the  size  of  production  terms  in  the  turbulence  transport  equations  upon  "cold-start" 
initialization. 

4.2  Turbulenl  Flow  Through  a  Supersonic  Compressor  Cascade  Operating  at  Unique 

Incidence  Condition 

The  first  cascade  to  be  investigated  is  the  PAV-1.5  supersonic  compressor  cascade 
tested  at  DFVLR  by  Schreiber  (1988).  This  pre-compression  blade  was  designed 
especially  to  investigate  shock-boundary  layer  interaction  with  separation.  At  the  test 
freestream  Mach  number,  a  standoff  leading  edge  shock  forms,  which  gives  rise  to  a 
separated  shock-boundary  layer  interaction  aft  of  mid  chord  on  the  suction  surface  of  the 


Figure  4.3  Boundary  layer  Reynolds  normal  stress  distributions  for  a  turbulent 
flat  pla^e  flow.  Experimental  data  due  to  Klebanoff  (1954) :  stream  wise 
intensity  (++++  symbols),  plate  normal  intensity  (oooo  symbols),  crossflow 
intensity  (xxxx  symbols).  Computations  using  hybrid  model  (solid  line). 
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adjacent  passage.  Though  the  measured  absolute  inlet  Mach  number  was  supersonic,  the 
blade  row  stagger  angle  was  high  so  the  axial  component  of  the  inlet  velocity  was 
subsonic.  This  gives  rise  to  the  "unique  incidence"  condition  wherein  there  exists  a  fixed 
relationship  between  inlet  Mach  number  and  inlet  flow  angle.  This  phenomena  as  well  as 
the  complex  wave  interaction  field  within  the  passage  and  shock-boundary-layer  interaction 
provide  a  challenging  test  case  for  both  numerical  scheme  and  turbulence  model. 

Prediction  of  the  unique  incidence  angle  requires  reasonably  good  resolution  of  both  the 
leading  edge  geometry  and  the  leading  edge  shock  structure. 

The  computed  case  was  experimentally  tested  in  air  at  an  inlet  Mach  number  of 
1.53  and  a  maximum  attainable  static  pressure  ratio  of  2.13.  The  measured  axial  velocity 
density  ratio  (AVDR)  of  1.02  indicated  that  the  flow  was  close  to  two-dimensional.  The 
Reynolds  number  based  on  chord  was  2.7  x  106.  The  inlet  turbulence  intensity  was 
measured  using  a  Laser-two-focus  (L2F)  velocimeter  to  be  no  more  than  1  %,  which  is  the 
value  used  in  the  computations.  As  mentioned  above,  the  inlet  Mach  number  is  supersonic, 
but  axial  velocity  at  the  inlet  to  the  computational  domain  is  subsonic  allowing  left  running 
characteristics  to  propagate  out  of  the  inlet  plane.  For  this  reason,  subsonic  inlet  boundary 
conditions  were  specified  :  pD  =  101325  N/m2,  T0  =  300  K,  V6oo=  379.5  m/s.  At  the 

subsonic  exit  plane  the  back  pressure,  pe  =  56500  N/m2,  was  specified  corresponding  to 
the  experimentally  measured  pressure  ratio  of  the  cascade,  p2/pi  =  2.13. 

The  129  x  100  computational  grid  used  was  generated  using  Sorenson’s  (1980) 
GRAPE  code,  modified  by  G^rski  (1983)  to  generate  H-grids,  and  is  shown  in  Figure  4.4, 
where  only  every  other  grid  line  in  the  £  and  r|  directions  are  plotted,  for  clarity.  The  blade 
normal  grid  spacing  at  the  wall  was  prescribed  as  .00001 1  chord.  This  yielded  values  of 
y*  =  1  at  grid  points  adjacent  to  the  walls.  Except  in  the  immediate  vicinity  of  the  leading 


L29  x  100  computational  grid  for  the  PAV-1.5  cascade. 
y  every  other  grid  line  is  shown  in  both  £  and  il  directions. 
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and  trailing  edges,  the  suction  and  pressure  surface  boundary  layers  had  at  least  9 
grid  points  with  values  of  y+  <  20. 

The  convergence  history  for  this  computation  is  shown  in  Figure  4.5.  It  took 
approximately  6500  iterations  for  this  calculation  to  converge  as  measured  by  the  invariance 
of  total  number  of  supersonic  gridpoints  in  the  field.  It  was  not  possible  to  "cold  start”  the 
initialized  flowfield  at  the  specified  pressure  ratio,  because  the  code  became  rapidly 
unstable  when  this  was  attempted.  Rather,  the  back  pressure  had  to  be  increased  in  a 
stepwise  fashion  with  iteration  (notice  "jumps”  in  convergence  at  iteration  500, 1000, 
2000),  until  the  experimentally  imposed  pressure  ratio  could  be  specified  at  iteration  2000. 
The  "unhealthy"  convergence  history  is  due  in  part  to  the  highly  clustered  grid  and  also  to 
the  nearly  choked  operating  condition. 

At  the  time  it  was  run,  the  solution  using  the  two-dimensional  version  of  the  code 
required  approximately  40  minutes  of  CPU  time  on  a  Y-MP.  For  this  cascade 
computation,  Aty  <  Atc  for  40.0%  of  the  grid  points  in  the  domain  at  convergence.  The 
artificial  dissipation  was  scaled  as  in  equation  [3.20],  where  the  reference  velocity,  Vref  = 
a^Moo,  was  used.  With  this  scaling,  a  significant  recirculation  zone  exists  aft  of  a  shock 
boundary  layer  interaction  region  aft  of  mid-chord  on  the  suction  surface.  This  region, 
extends  well  out  into  the  passage  (nj  =  25),  and  reattaches  a  short  distance  downstream. 
When  the  velocity  scaling  was  used  exactly  as  presented  in  Chapter  3,  significant  local 
oscillations  in  flowfield  properties  were  seen  to  develop  in  this  region,  immediately  prior  to 
the  code  blowing  up.  This  is  not  surprising  since  though  local  velocities  are  small  in  the 
recirculation  zone,  physical  diffusion  terms  in  the  momentum  equations  are  also  small. 
Therefore,  significant  artificial  dissipation  is  required  to  stabilize  the  solution  in  this  region. 
The  solution  was  stabilized  in  an  "ad-hoc"  fashion  for  this  calculation,  by  implementing  the 


velocity  scaling  given  in  equation  [3.20],  only  in  the  immediate  vicinity  of  the  blade 
surfaces  0  <  jspecified)- 

In  Figure  4.6,  a  hand  rendering  of  the  shock  wave  pattern  deduced  from  Schlieren 
visualization  has  been  reproduced  from  Schreiber  (1988),  alongside  the  computed  shock 
wave  pattern  presented  as  divergence  of  velocity  contours.  The  key  features  of  the 
flowfield  are  evident  in  this  diagram,  including  the  bow,  lamda  and  passage  shocks.  In 
both  experiment  and  computation,  the  bow  shock  is  seen  to  impinge  on  the  suction  surface 
boundary  layer  of  the  adjacent  passage.  This  gives  rise  to  a  lambda  shock  structure,  a  rapid 
thickening  and  separation  of  the  boundary  layer,  and  a  Mach  reflection  which  impinges  on 
the  pressure  surface  of  the  same  passage.  The  high  pressure  ratio  operating  condition  of 
this  test  case  gives  rise  to  a  normal  passage  shock  which  impinges  upstream  of  midchord 
on  the  pressure  surface.  This  feature  is  also  evident  in  both  experiment  and  computation. 
The  computation  also  shows  some  evidence  of  an  oblique  trailing  edge  shock,  typical  of 
supersonic  compressor  cascades  at  high  operating  pressure  ratios. 

In  Figure  4.7,  the  predicted  isentropic  blade  surface  Mach  number  is  plotted 
against  the  experimental  values.  The  calculation  and  experiment  show  fairly  good 
agreement  The  features  labelled  A,  B  and  C  in  Figure  4.7  correspond  to  local 
compression  regions  where  the  bow  shock  impinges  on  the  suction  surface,  the  Mach 
reflection  impinges  on  the  pressure  surface  and  the  passage  shock  impinges  on  the  pressure 
surface. 

In  Figure  4.8,  the  computed  total  pressure  ratio  is  compared  with  traverse  probe 
measurements  at  an  axial  location  0.09  chord  downstream  of  the  cascade  exit  plane.  The 
wake  profile  and  loss  distribution  is  reasonably  well  predicted,  with  the  losses  associated 
with  the  lambda  shock  system  underpredicted.  The  wake  centerline  total  pressure  ratio  is 


110 


Figure  4.6  Shock  wave  pattern  for  PA V  - 1.5  cascade.  Divergence  of 
velocity  contours  (-300  to  -4800  by  -500  [s1]).  In  this  diagram,  the  top  two 
passages  show  computed  contours.  The  bottom  passage  is  the  shock  wave 
pattern  deduced  from  flow  visualization  and  L2F  measurements,  adapted  from 

Schreiber  (1988). 
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Figure  4.7  Isentropic  blade  surface  Mach  numbers  and  shock  structure 
identification  for  PAV-1.5  cascade  computation.  Calculated  (solid  line)  and 
experimental  values  (symbols).  The  bottom  half  of  the  plot  shows  the  shock 
wave  pattern  deduced  from  flow  visualization  and  L2F  measurements,  adapted 

from  Schreiber  (1988). 
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predicted  reasonably  well  considering  the  difficulty  in  L2F  measurement  at  this  location.  It 
is  noted  that  the  results  presented  are  not  fully  grid  independent  Small  changes  in  the 
pitchwise  grid  clustering  near  midpassage  gave  rise  to  as  much  as  a  2%  chord  difference  in 
the  impingement  location  of  the  bow  shock  on  the  suction  surface  and  a  5%  chord 
difference  in  the  location  of  the  passage  shock.  The  loss  distribution  aft  of  the  blade  was 
hardly  affected,  but  the  blade  surface  Mach  number  distributions  varied  noticeably. 

In  Figure  4.9,  computed  velocity  and  turbulence  intensity  profiles  at  four  locations 
on  the  suction  and  pressure  surfaces  are  presented.  Along  the  suction  surface,  the 
predicted  boundary  layer  is  seen  to  remain  quite  thin  for  this  high  Reynolds  number  flow, 
thickening  to  about  .01  chord  at  midchord.  At  .75  chord  the  boundary  layer  has  separated 
due  to  the  bow  shock  impinging  at  .67  chord,  and  the  boundary  layer  thickness  is  seen  to 
have  rapidly  increased  to  approximately  .03  chord.  At  .90  chord,  the  separated  boundary 
layer  has  grown  to  .05  chord.  The  turbulence  intensity  profile  behaves  in  a  manner 
consistent  with  a  boundary  layer  separation.  Namely,  aft  of  the  onset  of  separation  the 
turbulent  kinetic  energy  boundary  layer  thickens  rapidly  and  the  peak  in  intensity  appears 
well  away  from  the  blade  surface.  Careful  examination  of  Figure  4.9  also  shows  that  the 
turbulence  intensity  has  been  amplified  well  outside  of  the  boundary  layer.  This 
amplification  is  presumably  due  to  the  influence  of  the  shock  on  the  normal  stress 
components  of  the  production  term  in  the  turbulent  kinetic  energy  equation. 

The  predicted  boundary  layers  along  the  pressure  surface  are  seen  to  remain  quite 
thin  along  the  entire  blade.  The  influence  of  the  passage  shock  at  .40  chord  separates  the 
boundary  layer  at  .50  chord,  but  the  flow  reattaches  and  the  boundary  layer  thickness 
remains  approximately  .02  chord  from  .75  to  .90  chord.  Similar  trends  are  noticed  in  the 
local  turbulence  intensity  profiles,  with  the  typical  peak  away  from  the  blade  just  aft  of 
separation,  returning  very  close  to  the  wall  some  distance  after  reattachment. 
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suction  surface 


pressure  surface 


Figure  4.9  Local  boundary  layer  velocity  and  turbulence  intensity  profiles  at 
four  chord  locations  along  the  suction  and  pressure  surfaces  for  the  PAV-1.5 
cascade  computation.  Locations  a,  b,  c,  d  correspond  to  .25,  .50,  .75,  .90 
chord  respectively.  The  hash  marks  correspond  to  the  estimated  boundary  layer 
thicknesses  reported  by  Schreiber  (1988). 


Schreiber  (1988)  provided  measured  loss  coefficients  at  maximum  attainable 
cascade  pressure  ratio  for  a  number  of  operating  inlet  Mach  numbers.  For  comparison,  the 
code  was  run  at  two  additional  operating  points  within  the  envelope  of  the  experimental 
tests  (P2/P1  =  2.01,  Miniet  =  1.45  and  P2/P1  =  2.28,  Minlel  =  1.58).  Figure  4.10  shows 
computed  total  pressure  loss  coefficients  at  all  three  operating  points  computed  along  with 
the  envelope  of  experimental  loss  coefficients.  Computed  values  lie  within  the  envelope  of 
experimental  values.  It  is  noted  that  Schreiber  attributes  the  scatter  in  measured  loss 
coefficient  to  variations  in  experimental  axial  velocity  density  ratio. 

Predicted  and  measured  performance  parameters  for  this  cascade,  operating  at  the 
design  condition  (P2/P1  =2.13,  M^j^  =  1.53)  are  presented  in  Table  4.1. 

Table  4.1.  Comparison  of  Performance  Parameters  for  PAV-1.5  Cascade 


Flowfield  Parameter 


Measured 

Computed 

Difference 

CLa 

=  0.43 

0.38 

=  12% 

G)b 

0.144 

.149 

3.4  % 

Pi 

-1.8° 

-1.3° 

0.5° 

P2 

2.7° 

3.1° 

0.4° 

S/C«PSS 

.63 

.67 

6.3  % 

S/C«Pps 

=  0.40 

.40 

=  0.0 

AVDR 

1.02 

1 

NA 

CL  = 


lift  coefficient  computed  from 


(Lift  per  unit  span)xym 
.5pmV2m 


Cl>  =  ; 


Po~  "  Po 


pressure  loss  coefficients  computed  from  Po~ '  P« 
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4.3  Turbulent  Flow  Through  a  Subsonic  Cascade  of  Double  Circular  Arc  Profiles  with 
Flow  Reversal 

The  second  cascade  flow  to  be  computed  is  the  Applied  Research  Laboratory 
(ARL)  double  circular  arc  cascade  tested  at  Penn  State  by  Zierke  and  Deutch  (1989).  The 
computed  case  was  tested  at  a  negative  incidence  of  1.5  degrees.  The  working  fluid  was 
air  at  standard  atmosphere  with  an  inlet  velocity  of  32.9  m/s  (inlet  Mach  number  =  0.1). 
The  Reynolds  number  based  on  chord  was  5.0  x  105.  Freestream  turbulence  intensity  was 
measured  at  1.8  %  at  20  %  chord  along  the  suction  surface,  which  is  the  value  used  in  the 
computation.  The  axial  velocity  ratio  was  measured  to  be  between  0.97  and  1.03, 
indicating  that  the  flow  was  close  to  two-dimensional. 

The  129  x  91  computational  grid  used  was  generated  using  the  GRAPE  code,  and 
is  shown  in  Figure  4.1 1.  Grid  spacing  in  the  blade  normal  direction  was  set  to  .000023 
chord  on  the  blade  surfaces.  This  yielded  values  of  y+  ~  1  at  grid  points  adjacent  to  the 
walls.  Except  in  the  immediate  vicinity  of  the  leading  and  trailing  edges,  the  suction  and 
pressure  surface  boundary  layers  had  at  least  1 1  grid  points  with  values  of  y+  <  20. 

It  was  only  possible  to  obtain  a  steady  solution  when  a  coarse  "preliminary"  grid 
was  used  for  this  case.  These  coarse  grid  calculations  overpredicted  skin  friction  along  the 
entire  length  of  the  suction  surface,  so  the  flow  remain  attached  and  a  steady  state  solution 
was  achieved.  The  more  refined  grid,  shown  in  Figure  4.1 1,  adequately  resolved  both 
inner  layer  and  core  flow  regions  yielding  more  accurate  skin  friction  and  boundary  layer 
profiles.  However,  because  both  calculation  and  experiment  show  regions  of  mean  flow 
reversal  near  the  trailing  edge,  it  was  not  possible  to  obtain  a  steady  solution. 
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The  convergence  history  for  this  computation  is  shown  in  Figure  4.12.  It  took 
approximately  7000  iterations  for  this  calculation  to  acquire  a  4.5  order  of  magnitude  drop 
in  the  RMS  density  residual.  At  the  time  this  computation  was  performed,  this 
corresponded  to  approximately  36  minutes  of  CPU  time  on  the  Y-MP.  However  as  shown 
in  Figure  4.12,  after  7000  iterations,  the  residual  changes  begin  to  increase  and  then  level 
off.  This  is  attributed  to  periodic  shedding  of  vorticity  from  the  aft  portion  of  the  suction 
surface.  For  this  cascade  computation,  A^  <  Atc  for  7.5  %  of  the  grid  points  in  the 
domain  when  the  calculation  was  stopped  at  14000  iterations. 

Despite  the  lack  of  a  steady  state  solution,  the  flow  along  the  blade  remained 
relatively  unchanged  after  7000  iterations  except  for  quasi-periodic  shifts  in  the  boundary 
layer  velocity  and  turbulence  intensity  profiles  near  the  aft  end  of  the  suction  surface.  The 
measured  flow  also  showed  a  small  region  of  mean  backflow  near  the  trailing  edge  of  the 
blade  (Zierke  and  Deutsch,  1989),  and  for  that  reason  was  also  probably  unsteady.  In 
Figure  4.13,  comparison  is  made  between  computed  blade  surface  pressure  coefficient  and 
measured  values.  Agreement  is  good  along  both  blade  surfaces.  The  oscillations  in  the 
pressure  distributions  near  the  leading  and  trailing  edges  in  Figure  4.13  are  caused  by  the 
velocity  scaling  of  the  artificial  dissipation.  The  artificial  dissipation  was  scaled  as  in 
equation  [3.20],  where  the  reference  velocity,  Vref  =  a*^.  The  H-grid  used  gives  rise  to 
highly  skewed  regions  near  the  relatively  blunt  leading  and  trailing  edges  of  this 
configuration,  causing  the  velocity  scaling  presented  "as  is"  to  give  rise  to  these 
oscillations. 

In  Figure  4.14,  the  predicted  boundary  layer  profiles  at  three  chordwise  locations 
on  the  suction  surface  are  plotted  with  those  measured  by  laser  Doppler  velocimeter. 
Agreement  is  excellent  at  20  %  chord  and  50%  chord  and  reasonable  at  90  %  chord.  Local 
turbulence  intensity  profiles  are  presented  for  three  chordwise  locations  on  the  suction 
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Figure  4.13  Pressure  coefficient  for  ARL  DCA  cascade  computation. 
Calculated  (solid  line)  and  experimental  values  (symbols). 


surface  in  Figure  4.15.  As  above,  agreement  between  calculation  and  experiment  is  good 
at  the  first  two  stations,  and  reasonable  in  the  aft  portion  of  the  blade. 


Predicted  and  measured  performance  parameters  for  this  cascade  are  presented  in 
Table  4.2. 


Table  4.2.  Comparison  of  Performance  Parameters  for  ARL-DCA  Cascade 


Flowfield  Parameter 

Measured 

Cl' 

0.82 

b 

CO 

0.094 

Pi 

-1.5° 

P2 

14.1° 

S/Csepss 

.45 

S^CscPps 

none 

AVDR 

0.97-1.03 

Computed 

Difference 

0.88 

7  % 

0.111 

18% 

prescribed 

NA 

13.0° 

1.1° 

fluctuates  near  0.90 

NA 
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NA 

1 

NA 

lift  coefficient  computed  from 
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pressure  loss  coefficients  computed  from 
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CHAPTER  5 

PREDICTION  OF  THREE-DIMENSIONAL  INTERNAL 
AND  TURBOMACHTNERY  FLOWS 

In  order  to  provide  three-dimensional  validation  of  the  turbulence  modelling  and 
numerical  methodologies,  the  code  was  used  to  predict  the  pressure  distribution,  primary 
velocity  and  secondary  flows  in  an  incompressible,  turbulent  curved  duct  flow  for  which 
CFD  validation  quality  data  is  available.  Also,  a  subsonic  compressor  rotor  passage,  for 
which  detailed  laser,  rotating  hot-wire  and  five-hole  pressure  probe  measurements  have 
been  made  is  computed.  Detailed  comparisons  between  predicted  and  measured  core  flow 
and  near  wall  velocity  profiles,  wake  profiles,  and  spanwise  mixing  effects  downstream  of 
the  rotor  passage  are  presented  for  this  case. 

5.1  Turbulent  Flow  in  a  90°  Bending  Duct  of  Square  Cross-Section 

The  ability  of  the  present  approach  to  predict  three-dimensional  turbulent  flows 
was  tested  by  computing  an  incompressible  duct  flow  investigated  experimentally  by 
Taylor  et  al.  (1981).  The  bulk  velocity  of  wafer  through  the  90°  bending  duct  of  square 
cross-section  was  1.00  m/s.  This  corresponds  to  a  Reynolds  number  based  on  hydraulic 
diameter  of  approximately  40000  and  a  Dean  number  of  approximately  18700.  A  short 
inlet  extension  to  the  duct  was  designed  by  the  experimenters  to  provide  thin  inlet  boundary 
layers  at  the  bend  entrance.  In  earlier  work  by  Humphrey  et  al.  (1981),  the  same  facility 
was  fitted  with  a  longer  inlet  section  so  as  to  provide  fully  developed  turbulent  velocity 
profiles  near  the  bend  entrance.  The  thinner  inlet  velocity  profiles  of  the  present  case  give 
rise  to  weaker,  but  significant,  secondary  motions  in  the  bending  portion  of  the  duct,  and 


therefore  this  case  provides  an  appropriate  three-dimensional  test  case  for  the  numerical 
scheme  and  the  turbulence  model. 

An  interesting  difficulty  in  computing  this  flow  arises  from  the  fact  that  the  velocity 
profiles  near  the  bend  inlet  are  significantly  influenced  by  the  development  of  secondary 
motions  of  the  second  kind  in  the  the  straight  duct  section  upstream  of  the  bend.  Accurate 
computation  of  secondary  flows  within  the  duct  relies  primarily  on  rotational  inviscid 
transport  of  accurately  prescribed  near  wall  inlet  transverse  vorticity  profiles  [see 
Lakshminarayana  and  Horlock  (1973)  for  example].  However,  isotropic  eddy  viscosity 
models  are  incapable  of  predicting  the  Reynolds  stress  anisotropies  which  give  rise  to 
secondary  motions  of  the  second  kind  and  therefore  to  the  "bowed"  inlet  velocity  profiles 
these  motions  give  rise  to.  This  difficulty  seems  to  discount  the  use  of  a  computational  grid 
with  an  inlet  extension  corresponding  to  the  experimental  apparatus,  if  the  present  two- 
equation  turbulence  model  is  used.  In  the  present  work,  this  difficulty  is  overcome  by 
setting  the  inlet  boundary  very  close  to  the  bend  inlet.  Measured  velocity  and  pressure 
distributions  are  used  to  provide  inlet  boundary  conditions  which  are  physically  plausible 
enough  to  allow  a  convergent  Navier-Stokes  solution  despite  the  proximity  of  the  inlet 
boundary  to  regions  where  significant  flowfield  gradients  exist. 

In  Figure  5.1,the81x41x21  computational  grid  used  for  this  calculation  is 
shown.  The  half  duct  was  computed  here  to  exploit  the  symmetry  of  the  problem.  Grid 
spacing  at  the  wall  was  0.00033  times  the  duct  height  which  gave  rise  to  wall  adjacent  grid 
point  values  of  y+  =  1.5  along  the  entire  length  of  the  duct.  The  inlet  extension  was  1/4 
hydraulic  diameter  upstream  of  the  bend  entrance.  Measured  primary  velocity  magnitudes 
and  the  radial  static  pressure  distribution  at  this  location  were  interpolated  onto  the  inlet  grid 
plane.  A  composite  law-of-the-wall  due  to  Spalding  [see  White  (1974)  for  instance]  was 
used  to  approximate  the  primary  velocities  between  the  wall  and  the  first  measurement 
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Figure  5.1  81  x  21  x  41  Computational  grid  and  notation  for 
incompressible  duct  flow  computation. 
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location  away  from  the  wall  (0.1  times  the  duct  height).  The  simulation  was  run  with  an 
average  inlet  Mach  number  close  to  0.2,  and  the  subsonic  Riemann  extrapolation  boundary 
condition  described  in  Chapter  3  was  used,  to  accommodate  the  compressible  formulation 
of  the  flow  solver.  Inlet  total  pressure  and  temperature  profiles  were  specified  from  the 
above  prescriptions  and  the  assumption  of  constant  inlet  static  temperature.  Additionally, 
zero  inlet  flow  angles  and  a  back  pressure  adjusted  to  pass  the  mass  flow  necessary  to 
achieve  Reynolds  number  similarity  were  specified.  Inlet  turbulence  intensity  was  set  to 
3%  corresponding  to  approximate  measured  inlet  average  intensity.  The  artificial 
dissipation  was  scaled  by  (Vi^/V^)2,  as  in  equation  [3.20], where  the  mean  inlet  velocity 
was  used  as  the  reference  velocity,  Vref. 

The  code  was  run  for  3000  iterations,  at  which  point  the  RMS  density  and 
turbulent  kinetic  energy  residuals  had  dropped  3  orders  of  magnitude.  The  computation 
was  then  stopped,  since  the  pressure  and  primary  and  secondary  velocity  fields  remained 
unchanged  after  about  2500  iterations.  The  convergence  history  for  this  computation  is 
shown  in  Figure  5.2.  This  computation  required  approximately  7  hours  on  an 
IBM  3090-600S.  At  convergence,  the  parabolic  stability  constraint  was  less  restrictive  than 
the  hyperbolic  constraint  for  all  of  the  grid  points  in  the  domain,  that  is  Aty  >  Atc 
everywhere.  This  is  due  to  the  fact  that  since  the  duct  walls  (and  hence  the  near  wall 
damping  effects  in  the  low-Reynolds  number  turbulence  model)  extend  from  inlet  to  exit, 
there  is  no  region  in  the  computed  flow  where  both  large  effective  diffusivities  and  large 
grid  metrics  occur.  Therefore,  the  first  term  in  equation  [3.46]  is  less  than  the  second  at 
every  grid  point. 

In  Figure  5.3,  the  predicted  streamwise  distribution  of  pressure  coefficient  along 
the  symmetry  plane  is  presented.  Agreement  between  the  experiment  and  the  computation 
at  all  five  radial  locations  is  very  good.  The  high  wave  number  pressure  oscillations  near 


Figure  5.3  Stream  wise  static  pressure  coefficient  distribution  at  several 
locations  along  symmetry  plane  for  incompressible  duct  flow  computation. 
Calculated  (solid  line)  and  experimental  values  (symbols). 


the  inlet  section  are  attributed  to  the  proximity  of  the  imposed  inlet  boundary  conditions  to 
the  bend. 
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Streamwise  velocity  contours  are  plotted  at  six  streamwise  locations  in  Figure  5.4, 
and  show  reasonably  good  agreement  with  the  experimental  data.  The  nearly  exact 
correspondence  at  the  inlet  (X/D  =  -0.25)  is  due  to  boundary  condition  specification  from 
the  data.  The  development  of  the  secondary  vortex  at  the  inner  wall  and  its  effect  on 
velocity  defect  in  the  mainstream  is  captured  accurately.  Both  the  size  of  this  vortex  and  the 
magnitude  of  the  streamwise  velocities  in  this  secondary  flow  region  are  well  predicted. 
The  secondary  velocities  are  shown  in  Figure  5.5.  Here  again,  agreement  is  good.  The 
ability  of  the  code  to  capture  the  far  downstream  flow,  where  considerable  redistribution  of 
the  flowfield  has  occurred  is  evident  from  both  Figures  5.4  and  5.5. 

The  good  agreement  between  experiment  and  computation  achieved  in  this 
validation  study  provides  evidence  that  the  numerical  and  modelling  strategies  used  can 
yield  acceptable  solutions  of  complex  three-dimensional  internal  flows. 


As  stated  in  Chapter  1,  the  ultimate  goal  of  this  thesis  work  is  to  provide  accurate, 
detailed  predictions  of  viscous  flows  in  complex  three-dimensional  turbomachinery 
configurations.  Accordingly,  the  code  has  been  utilized  to  predict  the  flow  field  in  a 
subsonic  rotor  operating  at  peak  pressure  rise  coefficient.  The  so  called  "Penn  State 
Compressor  Rotor"  has  a  hub  to  tip  ratio  of  0.5  with  an  outer  annulus  wall  diameter  of 
0.936  m.  The  rotor  has  21  blades.  The  chord  length  increases  from  0.124  m.  at  the  hub  to 
0.154  m.  at  the  tip.  The  blade  stagger  angle  increases  from  22.5°  at  the  hub  to  45.0°  at  the 
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Figure  5.4  Contours  of  normalized  streamwise  velocity,  V#/Vb,  for  incompressible  duct  flow 


computation.  Experimental  contours,  reproduced  from  Taylor  et  al.  (1981),  appear  on  bottom  side  of  each 
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Figure  5.5  Secondary  velocities  at  various  streamwise  stations  for  incompressible  duct  flow 
computation,  a)  Computed  secondary  velocity  vectors  [solid  lines  indicate  location 
of  traverses  in  b)].  b)  Scaled  radial  velocities,  V/Vb,  at  five  radial  locations. 
Calculated  (solid  line)  and  experimental  values  (symbols).  See  Figure  5.1  for  notation. 
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tip.  The  tip  clearance  was  either  2  mm  or  3.3  mm  for  the  various  experimental  tests  as 
specified  below. 

The  operating  conditions  for  the  computed  case  are  as  follows  :  flow  coefficient, 

O  =  0.50,  pressure  rise  coefficient  (peak),  V  =  0.55,  rotor  shaft  angular  velocity,  co  =  1 13 
S'1,  tip  velocity,  U  =  53  m/s.  The  Reynolds  number  of  the  flow  based  on  chord  length  at 
midspan  and  rotor  tip  speed  was  480,000. 

A  large  amount  of  laser  Doppler  velocimeter,  hotwire  and  five  hole  pressure  probe 
data  is  available  for  this  rotor.  The  axial  and  tangential  components  of  relative  velocity 
were  measured  from  blade  to  blade  at  various  chordwise  locations  inside  the  passage,  using 
an  ensemble  averaging  technique  (Popovski  and  Lakshminarayana  (1986)).  The  flowfield 
near  the  blade  surfaces  was  measured  using  a  two  sensor  miniature  hotwire  probe,  rotating 
with  the  rotor  (Lakshminarayana  and  Popovski  (1987)).  Additionally,  the  wakes 
downstream  were  measured  recently  by  Prato  (1990)  using  a  five  hole  pressure  prc’oe.  For 
the  passage  flow  measurements,  the  tip  clearance  was  2  mm.  The  blade  tips  were  shaved 
to  provide  a  3.3  mm  tip  clearance  prior  to  the  wake  measurements. 

In  Figure  5.6a,  the  89  x  45  x  45  (180225  points)  computational  grid  used  for  this 
calculation  is  shown.  This  H-grid  was  generated  using  a  interactive  program  developed  by 
Beach  (1990),  which  is  described  in  detail  there.  The  normal  grid  spacing  on  the  blade 
surfaces  was  4.7  x  10"4  span.  The  normal  grid  spacing  at  the  hub  and  shroud  was  1.9  x 
10'3  span.  This  yielded  values  of  y+=  6  at  grid  points  adjacent  to  the  blade  surfaces  and 
y+=  30  along  the  endwalls.  It  is  noted,  that  for  accurate  implementation  of  nearwall  two- 
equation  turbulence  models,  values  of  y+  should  be  closer  to  1  at  the  first  gridpoint  off  the 
wall.  Resource  limitations  allowed  only  45  blade  to  blade  and  hub  to  casing  grid  points. 
For  this  high  Reynolds  number  rotor  flow,  the  resulting  grid  generation  compromised 
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Figure  5.6  Computational  grid  and  configuration  details  for  PSU  rotor  flow  computation,  a)  89  x  45 


45  computational  grid,  b)  detail  of  "thin  blade"  tip  clearance  gap  at  midchord,  c)  hub-leading  edge 


detail  showing  relative  flow  velocities  at  stationary-rotating  hub  interface  in  relative  frame. 
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between  refined  nearwall  and  core  region  grid  resolution  in  order  to  provide  reasonable  grid 
clustering  everywhere. 

In  order  to  account  for  tip  clearance  flow,  the  grid  was  modified  in  the  immediate 
vicinity  of  the  shroud,  as  shown  in  Figure  5.6b.  This  allows  for  three  periodic  flow 
through  grid  points  in  the  gap  region.  Such  a  "thin  blade"  tip  clearance  approach  has  been 
undertaken  by  several  researchers  [see  Bansod  and  Rhie  (1990)  for  example].  In  the 
compressor  facility,  the  location  of  the  interface  between  rotating  and  stationary  segments 
of  the  hub  is  3  mm  upstream  of  the  leading  edge  at  the  root.  To  accommodate  this  in  the 
flow  computation,  a  no-slip  relative  velocity  of  -corhub  is  enforced  at  all  grid  points  on  the 
hub  for  which  x  <  Xjje  -  3  mm.  This  is  illustrated  in  Figure  5.6c.  Measured  inlet  axial  and 
tangential  velocities,  provided  in  Popovski  and  Lakshminarayana  (1986),  were  interpolated 
onto  the  inlet  plane  to  provide  the  rotor  flow  inlet  boundary  conditions  discussed 
previously.  A  composite  law-of-the-wall  due  to  Spalding  [see  White  (1974)  for  instance] 
was  used  to  approximate  the  total  velocities  between  the  endwalls  and  the  first  measurement 
location  away  from  these  walls.  The  inlet  turbulence  intensity  was  specified  as  .02,  and  the 
inlet  turbulence  length  scale  as  .01  times  the  blade  spacing  at  midspan. 

The  code  was  run  for  5200  iterations,  at  which  point  the  RMS  density  and 
turbulent  kinetic  energy  residuals  had  dropped  5.0  orders  of  magnitude.  The  computation 
was  then  stopped,  since  the  streamwise  and  radial  near-wall  velocity  profiles  remained 
unchanged  after  about  4500  iterations.  The  convergence  history  for  this  computation  is 
shown  in  Figure  5.7.  This  computation  required  approximately  25  hours  on  a  Cray  2.  At 
convergence,  the  parabolic  stability  constraint  was  more  restrictive  than  the  hyperbolic 
constraint,  that  is  Aty  <  Atc,  for  6.5  %  of  the  grid  points  in  the  domain.  The  artificial 
dissipation  was  scaled  according  to  equation  [3.20],  where  the  rotor  tip  speed,  U,  was 
used  as  the  reference  velocity. 
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Passage  Flowfield 

In  Figure  5.8,  the  blade  to  blade  axial  velocity  distributions,  obtained  at  several 
spanwise  locations  by  Popovski  and  Lakshminarayana  (1986),  using  LDV  measurements, 
is  presented  with  the  Navier-Stokes  predictions.  The  rotational  inviscid  flow  features  are 
seen  to  be  well  captured  by  the  simulation,  including  immediately  upstream  of  the  leading 
edge  where  the  presence  of  the  blade  gives  rise  to  small  velocity  magnitudes  due  to  the 
proximity  of  the  stagnation  point  (inviscid  effect).  The  effect  of  tip  clearance  flow,  which 
interacts  with  the  core  flow  near  the  tip  (r/r^p  =  0.96)  to  produce  a  velocity  defect  near  the 
midpassage  is  captured  qualitatively  by  the  code. 

Detailed  blade  boundary  layer  measurements  were  made  for  this  operating 
condition,  using  a  rotating  hot  wire,  by  Lakshminarayana  and  Popovski  (1987).  In 
Figures  5.9  and  5.10,  the  predicted  streamwise  and  radial  velocity  profiles  are  compared 
with  measured  values.  The  streamwise  boundary  layer  velocity  profile  near  the  hub  and 
midspan  regions  is  captured  with  reasonable  accuracy.  The  region  r/rlip  =  0.92  -  0.96 
shows  significant  discrepancy,  in  the  near  wall  region.  The  measured  profiles  in  the  near 
wall  region  of  the  tip  sections,  r/r^p  =  0.92  -  0.96,  show  a  tendency  to  separate,  and  this  is 
not  captured  well  by  the  simulation.  This  is  a  region  of  large  radial  boundary  layer 
transport  interaction  with  the  leakage  and  secondary  flow.  Considering  these  complexities, 
the  predictions  are  reasonably  good. 

The  radial  velocity  profiles  are  shown  in  Figure  5.10.  The  agreement  between 
experiment  and  computation  is  good  up  to  r/r^  =  0.92,  but  deteriorates  as  the  casing  is 
approached.  Some  of  the  discrepancy  is  attributed  to  poor  resolution  of  tip  clearance  flow 
and  the  resulting  blade  unloading  in  this  region.  The  radial  velocities  are  strongly 
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Figure  5.9  Suction  surface  stream  wise  boundary  layer  velocity  distribution 
(Wj/Ws  edge)  for  PSU  rotor  flow,  at  four  spanwise  locations.  Rotating  hot  wire 
measurements  due  to  Laxshminarayana  and  Popovski  (1987)  (symbols)  and 
computed  distributions  (solid  lines). 


Figure  5.10  Suction  surface  radial  boundary  layer  velocity  distribution  (W,/Ws 
edge)  for  PSU  rotor  flow,  at  four  spanwise  locations.  Rotating  hot  wire 
measurements  due  to  Lakshminarayana  and  Popovski  (1987)  (symbols)  and 
computed  distributions  (solid  lines). 
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dependent  on  both  the  streamwise  velocity  defect  near  the  blade  surface,  and  the  local  radial 
pressure  gradient.  Good  prediction  of  radial  velocities  at  r/r ^  =  0.58  and  0.75  is  due  to 
good  prediction  of  streamwise  velocity  profiles  in  this  region  (compare  Figure  5.9).  The 
discrepancy  between  radial  velocity  data  and  prediction  from  r/rd  =  0.92  to  0.96  is  due  to 

aw , 

poor  prediction  of  streamwise  velocity.  The  predicted  nearwall  gradients,  (Figure 

5.9) ,  in  these  regions  are  smaller  than  their  measured  values.  This  gives  rise  to  large 
discrepancies  between  measured  and  predicted  radial  velocities  in  the  tip  regions  (Figure 

5.10) . 

Rotor  Wakes 

Recently,  Prato  (1990)  obtained  extensive  wake  flow  measurements  for  this  rotor 
flow.  In  Figures  5.1 1  to  5.13,  measured  axial,  relative  tangential  and  radial  velocity 
profiles  in  the  near  and  far  wake  are  compared  with  computed  profiles.  Agreement 
between  calculated  and  measured  profiles  is  very  good.  Comparisons  are  included  only  for 
three  radius  ratios,  r/r^  =  0.57, 0.73  and  0.86.  Agreement  in  the  near  casing  region,  r/rtip 
=  0.96,  was  not  good,  due  to  the  dominant  tip  vortex  physics  in  this  region  and  the 
relatively  crude  gap  modelling  undertaken  herein.  The  simulated  wake  thicknesses  and 
spreading  rates  are  good.  Wake  centerline  velocities  are  all  well  predicted,  as  is  wake 
spreading  rate.  The  fairly  accurate  prediction  of  wake  radial  velocity  profiles,  as  seen  in 
Figure  5.13,  is  encouraging,  considering  the  low  relative  magnitude  of  these  spanwise 
flows.  The  wake  data  at  r/r^p  =  0.57  shows  the  influence  of  secondary  flow,  with  both 
inward  and  outward  radial  velocities.  Blade  rotation  gives  rise  to  radially  outward 
components  along  both  the  suction  and  pressure  surface  boundary  layers,  while  near  the 
hub,  spanwise  secondary  flow  components  are  radially  outward  on  the  suction  side  and 
inward  on  the  pressure  side.  Thus  the  two  effects  are  additive  on  the  suction  side  and 
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Figure  5.1 1  Wake  axial  velocity  distribution  (Wx/Wx  ^j^)  for  PSU  rotor  flow,  at  three  spanwise  locations. 
Five  hole  probe  measurements  due  to  Prato  (1990)  (symbols)  and  computed  distributions  (solid  lines). 
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Figure  5.12  Wake  tangential  velocity  distribution  (W^AVq  for  PSU  rotor  flow,  at  three 
span  wise  locations.  Five  hole  probe  measurements  due  to  Prato  (1990)  (symbols)  and  computed 

distributions  (solid  lines). 
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r/nip  =  0.57 


r/nip  =  0.73 


opposing  on  the  pressure  side.  This  effect  has  been  measured  as  well  as  predicted  as  seen 
in  the  wake  profiles  presented  in  Figure  5.13  (r/rup  =  0.57,  s/c=1.02). 


Spanwise  Mixing  Effects  and  Losses 

In  Figure  5.14,  spanwise  distributions  of  circumferentially  averaged  axial  velocity 
and  absolute  flow  angle  are  presented  at  five  downstream  locations,  ranging  from  the  very 
near  wake  to  the  far  wake  regions.  Agreement  between  experiment  and  computation  is 
excellent  except  near  the  tip.  The  simulation  captures  well  the  radial  distribution  of  these 
quantities  through  the  extent  of  the  wake.  Also,  the  streamwise  increase  in  Vx  and  decrease 
in  absolute  flow  angle  near  both  endwalls  is  captured. 

In  Figure  5.15,  the  computed  contours  of  relative  stagnation  pressure  loss 
coefficient  are  presented  for  grid  planes  coincident  with  the  blade  trailing  edge  and  the  exit 
of  the  computational  domain  in  the  far  wake.  A  number  of  sensible  and  interesting  flow 
features  are  apparent  in  these  plots.  The  loss  contours  at  the  trailing  edge  plane  (Figure 
5.15a)  reveal  the  presence  of  boundary  layers  near  the  blade  surfaces  where  peak  loss 
coefficients  as  high  as  0.50  are  predicted.  The  presence  of  secondary  flow  in  the  hub  wall 
region  and  the  resulting  migration  of  low  momentum  fluid  toward  the  suction  side  is 
indicated  by  the  presence  of  a  large  loss  region  near  the  intersection  of  the  suction  surface 
and  hub  wall.  The  presence  of  high  losses  near  the  outer  wall  reveal  the  qualitative  effect 
of  leakage  flow.  The  scraping  vortex  (of  sma'l  radius)  at  the  intersection  of  the  pressure 
surface  and  the  annulus  wall  is  also  captured  qualitatively  by  the  computation. 

The  loss  contours  far  downstream  (Figure  5. 15b)  reveal  an  interesting  phenomena 
related  to  spanwise  and  lateral  mixing  downstream  of  the  blade  row.  The  wake  spreading 
(compare  Figures  5.1 1  to  5.13),  decreases  losses  near  the  wake  centerline  and  increases 
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Figure  5. 14  Passage  averaged  (mass  weighted)  values  of  a)  axial  velocity  and 
b)  flow  angle  downstream  of  PSU  rotor  passage.  Averaged  five  hole  probe 
measurements  due  to  Prato  (1990)  (symbols)  and  computed  distributions  (solid 
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losses  away  from  it.  The  effect  of  transverse  mixing  increases  losses  substantially 
everywhere  except  in  the  inviscid  core  region.  The  inviscid  core  region  is  located 
approximately  in  the  middle  third  of  the  passage,  where  the  losses  are  nearly  zero.  The 
comer  flow  loss  region  (at  the  intersection  of  suction  surface  and  hub  wall)  has  grown  both 
in  the  blade  to  blade  and  spanwise  direction.  The  radial  velocities  in  the  wake  on  the 
suction  side  have  transported  and  elongated  this  loss  region  to  higher  radius. 

An  interesting  phenomenon  is  identified  in  the  tip  region  of  Figure  5.15b.  Here 
the  losses  have  not  only  decreased,  but  very  close  to  the  wall  negative  loss  coefficients  are 
predicted.  This  is  caused  by  spanwise  and  lateral  mixing  due  to  wakes,  leakage  vortices 
and  secondary  flows.  This  is  consistent  with  the  spanwise  mixing  analysis  carried  out  by 
Adkins  and  Smith  (1982),  whose  results  have  been  reproduced  in  Figure  5.16.  Their  input 
loss  distribution  was  based  on  empirical  correlations.  The  computation  was  based  on  a 
spanwise  mixing  analysis.  The  measured  values  indicate  negative  loss  coefficient  near  the 
outer  wall,  consistent  with  the  present  computational  results.  Additionally,  the  PSU  rotor 
computation  predicts  that  the  loss  core  in  the  annulus  wall  region  is  washed  out  by  mixing. 
This  is  clear  evidence  that  the  high  loss  region  observed  at  the  trailing  edge  in  the  annulus 
wall  region  is  "peeled  off'  by  secondary  flows  and  the  entire  loss  core  region  is  convected 
and  diffused  to  inner  radii.  The  spanwise  mixing  effects  are  thus  simulated  by  this  Navier- 
Stokes  code,  supporting  the  physical  phenomena  postulated  by  Adkins  and  Smith. 

The  passage  averaged  (mass  weighted)  losses  just  downstream  of  the  trailing  edge, 
both  measured  by  Prato  (1990)  and  predicted  by  the  present  authors,  are  shown  in  Figure 
5.17a.  The  profile  losses,  away  from  the  endwalls  are  captured  with  reasonable  accuracy. 
Underprediction  of  losses  at  hub  and  midspan  regions  are  consistent  with  the 
underprediction  of  wall  shear  stresses  along  the  blade  surfaces,  evident  upon  close 
examination  of  Figure  5.9.  The  data  includes  one  point  in  the  tip  region  and  the  prediction 
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Figure  5.16  Loss  coefficient  distributions  for  General  Electric  Low  Speed 
Research  Compressor.  This  figure  is  reproduced  from  Adkins  and  Smith 
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Figure  5.17  Passage  averaged  (mass  weighted)  relative  stagnation  pressure  loss 
coefficient,  Cp  lo5S,  a)  in  the  near  wake  and  b)  .66  chord  downstream  of  hub  trailirg 
edge.  Averaged  five  hole  probe  measurements  due  to  Prato  (1990)  (symbols)  and 

computed  distributions  (solid  lines). 
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is  poor  there.  The  losses  predicted  approximately  one  chordlength  downstream  are  shown 
in  Figure  5.17b.  The  span  wise  mixing  phenomena,  referred  to  above,  decrease  losses  in 
the  tip  region,  and  this  feature  is  captured  by  the  code.  Major  redistribution  of  losses  occur 
only  in  the  tip  region,  where  intense  spanwise  and  lateral  mixing  due  to  leakage  flow, 
secondary  flow  and  wake  radial  velocities  occur.  This  is  consistent  with  the  results  of 
Adkins  and  Smith  (1982). 

Overall  Perfoimance 

Computed  and  measured  flow  coefficients  and  pressure  rise  coefficients  are 
presented  in  Table  5.1. 
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Table  5.1.  Comparison  of  Performance  Parameters  for 
Penn  State  Compressor  Rotor 


Flowfield 

Parameter® 

Measured*3 

Computed 

% 

Difference 

Flow 

Coefficient, 

<t> 

0.50 

0.495 

1.0C 

Pressure 

Rise 

Coefficient, 

V 

0.55 

0.563 

2.4 

_ pvxp0  dA 

-  fVl  l^  u^1  ’Po"“7 

<j)=J^LdA/A,  2^  j  pvxdA 

b  Values  from  Lakshminarayana  and  Popov  ski  (1987) 

c  Note  that  close  agreement  for  ^  is  due  to  choice  of  inlet  boundary  condition  (see 
discussion  above). 

It  it  evident  that  the  overall  pressure  rise  coefficient  (Table  5.1)  as  well  as  overall 


losses  (Figure  5.17)  are  reasonably  well  predicted. 
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CHAPTER  6 

TRANSONIC  CENTRIFUGAL  COMPRESSOR  COMPUTATIONS 

In  this  chapter,  results  of  computations  of  the  flowfield  within  a  modem  transonic 
centrifugal  compressor  stage  are  presented  and  interpreted.  These  results  are  compared 
with  available  L2F  meridional  velocity  measurements  within  the  impeller  passage,  shroud 
static  pressure  measurements  and  overall  performance  measurements.  First,  solutions 
using  the  low  Reynolds  number  k-E  model  are  provided.  The  hybrid  k-e  /  ARS  model  was 
also  used  to  compute  this  flow,  both  with  and  without  incorporation  of  the  production  by 
system  rotation  terms  (Rjj  in  equation  [2.14])  included  in  the  model,  and  results  of  these 
computations  are  also  included. 

6.1  Background 

The  centrifugal  compressor  chosen  for  detailed  computational  study  was  designed 
and  tested  by  Krain  et  al.  (1988, 1989).  The  low  specific  speed  impeller  (Ns  =  80)  has  24 
straight  blade  surfaces  (hub  to  tip)  to  allow  flank  milling  manufacture,  and  a  30° 
backsweep.  In  Figure  6.1,  a  photograph  of  the  impeller  is  reproduced  from  Krain  (1988). 
In  Figure  6.2,  the  periphery  of  the  computational  grid  used  in  the  simulations  is  provided. 
Except  for  the  modem  impeller  wheel,  the  experimental  facility  used  by  Krain  was  the  same 
as  that  used  by  Eckardt  (1975, 1976),  with  a  constant  area  vaneless  diffuser  to  allow  a 
wide  operating  range.  As  seen  in  Figure  6.2,  the  diffuser  flowfield  was  also  computed.  In 
Figure  6.3,  experimentally  obtained  performance  maps  for  the  stage  have  been  reproduced 
from  Krain  (1988),  and  the  operating  point  computed  is  indicated.  The  design  point  was 
chosen  for  the  present  calculations.  The  flow  parameters  corresponding  to  this  operating 
condition  were  :  mass  flow  rate  =  4.0  kg/s  (.167  kg/s  each  passage),  maximum  inlet 
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Figure  6.1  Impeller  wheel  tested  experimentally  by  Krain  (1988),  and  investigated 
computationally  in  this  chapter.  This  figure  was  reproduced  from  Krain  (1988). 
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Figure  6.3  Performance  maps  for  Krain  impeller,  indicating  computed  operating  point. 
This  figure  was  adapted  from  Krain  (1988). 
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relative  Mach  number  =  0.8  (at  shroud),  machine  rotation  rate  =  22,363  rpm, 

P0, diffuser  exit/PO, impeller  inlet  -  ^  • 

In  Figures  6.4  and  6.5,  meridional  and  streamwise  views  of  the  59  x  27  x  27 
computational  grid  (4301 1  grid  points)  used  for  this  calculation  are  shown.  Also 
appearing  in  Figure  6.4  are  the  locations  of  the  experimental  laser  planes  used  by  Krain  and 
referred  to  in  the  results  to  follow.  The  experimentally  measured  impeller  wheel  was  fitted 
with  a  spherical  nose  cone.  For  the  computations,  a  cusped  spinner  cone  was 
implemented,  as  shown  in  Figure  6.4,  to  accommodate  the  H-grid  topology  of  the  code. 

One  of  the  few  practicalities  in  computing  centrifugal  compressor  flowfields  arises  from  the 
fact  that  the  impeller  flowfield  is  dominated  by  shear  layers  which  encompass  the  entire 
passage.  For  this  reason,  nearly  grid  independent  solutions  were  obtained  in  the  present 
study,  by  using  27  grid  points  in  the  blade  to  blade  and  hub  to  tip  directions.  This 
phenomenon  has  also  been  noted  by  Hah  and  Krain  (1989)  (28704  grid  points),  and  is 
apparently  the  reason  for  the  relatively  coarse  cross  stream  grids  used  in  centrifugal 
compressor  calculations  by  Rhie  et  al.  (1985)  (14250  grid  points),  Dawes  (1988)  (17051 
grid  points)  and  Choi  and  Knight  (1989)  (81600  grid  points). 

The  grid  was  generated  algebraically  using  the  impeller  geometry  provided  by 
Krain  and  Hoffmann  (1989, 1990),  and  a  grid  clustering  routine  provided  by  Basson 
(1991).  Krain  (1990)  measured  the  impeller  inlet  and  exit  tip  clearance  at  design  operating 
speed  to  be  0.5  mm  and  0.2  mm,  respectively.  The  tip  gap  was  varied  linearly  along  the 
shroud  from  inlet  to  exit  in  the  grid  generation  procedure.  As  with  the  rotor  flow 
computation  presented  in  Chapter  5,  a  "thin  blade"  tip  gap  approximation  is  used  in  the 
present  calculations.  As  shown  in  Figure  6.5,  eight  grid  points  were  used  to  resolve  the 
gap  region.  Additionally,  the  leading  and  trailing  edges  were  cusped  over  two  streamwise 
grid  points.  Wall  normal  grid  spacing  was  4.4  x  I0-4  times  the  inlet  span  along  the  hub 
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Diffuser  Exit  (i  =  59) 


Impeller  Exit  / 
Diffuser  Inlet  (i  =  49) 


Figure  6.4  Meridional  view  of  59  x  27  x  27  computational  grid  for  Krain  impeller 
computation.  Streamwise  grid  indexing  and  location  of  experimental  laser  planes  are 

indicated  for  reference. 


Figure  6.5  Cross  stream  views  of  59  x  27  x  27  computational  grid  for  Krain  impeller 
computation,  a)  Impeller  inlet  grid  slice,  b)  grid  slice  near  impeller  exit  and  c)  tip  gap 

detail  at  inlet  grid  plane. 
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and  2.2  x  10-4  times  the  inlet  span  on  the  blade  surfaces.  This  clustering  yielded  converged 
values  of  y+  ranging  from  2  to  7  along  the  hub  endwall  and  both  biad  surfaces  at 
midchord,  for  the  k-e  solution  provided  below. 

As  discussed  previously,  curvature,  rotation  and  viscous  physical  phenomena 
strongly  influence  the  impeller  fowfield,  especially  in  the  aft  portions  of  the  passage. 
Specifically,  strong  secondary  motions,  caused  by  turning  of  blade  and  endwall  boundary 
layers,  and  coriolis  forces,  influence  the  flow  development  as  illustrated  schematically  in 
Figure  6.6a.  In  this  figure,  the  secondary  flow  effects  labelled  A  through  D  correspond  to  : 

A)  Secondary  flows  arising  from  axial  to  radial  turning  of  blade  boundary  layers. 

B)  Secondary  flows  arising  from  axial  to  tangential  turning  of  endwall  boundary  layers. 

No  direction  convention  is  provided  for  this  effect,  since  endwall  boundary  layers  are 
convected  towards  the  suction  side  due  to  tangential  to  axial  turning  in  the  inducer,  but  axial 
to  tangential  turning  due  to  backsweep  of  the  impeller  under  consideration  gives  rise  to  the 
opposite  effect  in  aft  portions  of  the  passage. 

Figure  6.6b  is  a  sketch  which  illustrates  the  nature  of  these  curvature  induced  secondary 
motions.  Specifically,  the  three  dimensional  cross  flow  associated  with  axial  to  radial 
turning  of  the  pressure  side  blade  boundary  layer  (A  above)  is  illustrated. 

— •  — ♦ 

C)  Coriolis  "pipe"  force  or  "eddy"  arising  from  "2^  X  We  component  of  coriolis 
acceleration  (gives  rise  to  slip  velocity;  see  performance  discussion  below). 

D)  Coriolis  "spin"  force  arising  from  -2co  X  Wr  component  of  coriolis  acceleration. 

The  dominant  flow  feature  in  centrifugal  compressors  is  the  so  c^  d  "jet-wake" 
structure  which  forms  near  the  shroud  in  all  impeller  flows.  Cross  flows  which  are 
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Figure  6.6  Qualitative  representation  of  secondary  motions  present  in  a  centrifugal 
compressor  impeller. 

a)  Sense  of  curvature  and  rotation  induced  motions  : 

A)  Secondary  flows  arising  from  axial  to  radial  turning  of  blade  boundary  layers. 

B)  Secondary  flows  arising  from  axial  to  tangential  turning  of  endwall  boundary  layers. 
Q  Coriolis  "pipe"  force. 

D)  Coriolis  "spin"  force. 

b)  qualitative  sketch  of  pressure  surface  crossflow  arising  from  axial  to  radial  turning  of  the 
meridional  flow. 
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induced  by  axial  to  radial  turning  of  the  blade  boundary  layers  (effect  A  above),  accumulate 
low  momentum  boundary  layer  fluid  into  the  shroud  region  near  mid-passage.  Typically, 
the  pressure  side  vortex  is  stronger  than  the  suction  side  counterpart,  of  this  vortex  pair, 
labelled  A  in  Figure  6.6a.  For  this  reason,  the  location  of  this  region  of  low  energy  fluid  is 
typically  closer  to  the  suction  surface.  As  discussed  by  Eckardt  (1976)  and  others,  this 
pocket  of  low  momentum  fluid  forms  approximately  halfway  through  the  impeller  channel. 
It  typically  grows  to  fill  a  significant  portion  of  the  midpassage  region,  though  accounts  for 
a  significantly  smaller  percentage  of  total  mass  flow.  This  "wake"  flow  introduces 
significant  blockage  and  thereby  acceleration  of  the  surrounding  high  momentum  fluid, 
termed  the  "jet"  flow.  In  the  impeller  investigated  experimentally  by  Eckardt,  this  jet-wake 
structure  persisted  through  to  the  exit  of  the  impeller. 

Additionally,  tip  clearance  flow  can  significantly  alter  the  flow  pattern  in 
centrifugal  compressor  impellers.  Specifically,  some  10  %  of  the  total  mass  flow  passes 
through  the  tip  gap,  driven  by  the  pressure  gradient  between  suction  and  pressure  surfaces. 
Typically  this  leakage  flow  gives  rise  to  a  strong  vortex,  whose  sense  is  opposite  the  main 
flow  direction,  and  which  strays  into  the  adjacent  passage  and  interacts  with  the  wake  flow. 

Lastly,  the  direct  influence  of  streamline  curvature  and  system  rotation  on  the 
structure  of  local  turbulence  has  been  attributed  to  the  formation  and  persistence  of  the 
wake.  Specifically,  concave  curvature  along  the  shroud  tends  to  stabilize  the  turbulence 
locally,  which  should  in  turn  lessen  the  resistance  of  the  developing  shroud  boundary  layer 
fluid  from  being  transported  away  from  the  shroud  towards  the  hub  (transverse  flow 
separation).  Additionally,  the  relatively  distinct  boundary  often  observed  between  wake 
and  jet  flow,  may  be  attributable  to  locally  reduced  turbulence  mixing.  This  reduction 
could  arise  from  damping  near  the  shroud  due  to  concave  curvature  and/or  damping  near 
the  suction  surface  due  to  rotation  in  the  radial  portion  of  the  impeller. 
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In  short,  the  variety  of  physics  which  begin  to  manifest  themselves  in  the  aft 
sections  of  an  impeller  are  very  complex.  Since  the  performance  of  the  machine  is  critically 
dependent  on  these  phenomena,  especially  the  nature  of  exit  flow,  full  Navier-Stokes  and 
turbulence  modelling  methods,  which  can  provide  more  accurate  prediction  and 
understanding  of  these  phenomena,  are  useful  to  the  engineer. 


Often  in  unswept  impeller  computations,  relative  velocity  vectors  are  used  to 
illustrate  the  development  of  secondary  motions.  In  machines  with  blade  sweep,  the  choice 
of  an  appropriate  viewing  projection  can  be  difficult,  and  will  not  in  general  align  with  a 
local  grid  slice.  Accordingly,  velocity  vector  plots  are  not  very  useful  in  interpreting  the 
secondary  motions  in  the  impeller  under  consideration  here.  Normalized  helicity  density  is 
used  in  this  chapter  to  help  interrogate  the  cross  stream  physics. 


H  = 


W  c 


Normalized  relative  helicity,  lwl  M ,  is  a  useful  parameter  in  three- 
dimensional  postprocessing  of  vortical  flow  solutions.  Levy  et  al.  (1990)  recently  itemized 
the  advantages  of  normalized  helicity  as : 

1)  Vortices  can  be  distinguished  from  boundary  layers  (which  have  large  vorticity  but 
W  •  C  is  small). 


2)  Helicity  is  a  scalar,  and  as  such  does  not  introduce  the  ambiguities  mentioned  above  that 
velocity  and  vorticity  vectors  give  rise  to. 

3)  The  sign  of  the  helicity  indicates  the  rotational  sense  of  a  vortex. 

4)  Maxima  in  normalized  helicity  correspond  to  vortex  cores,  and  as  such  provide  a 
convenient  means  for  tracking  developing  vortices. 
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6.2  Flowfield  Calculation  Using  the  k-E  Model 

The  first  set  of  results  to  be  presented  are  for  a  solution  obtained  using  the  low 
Reynolds  number  k-e  model.  The  standard  quasi-one-dimensional  flow  initialization 
procedure  is  a  very  poor  initial  guess  for  this  flowfield,  since  centrifugal  pressure  rise 
dominates  the  rotational  inviscid  physics  in  this  machine,  but  is  not  accounted  for  in  such 
an  initialization.  Only  by  running  the  half  grid  problem  (30  x  14  x  14)  and  enforcing  inlet 
mass  flow  rate,  through  the  course  of  iteration,  was  it  possible  to  obtain  a  converged 
solution  from  the  quasi-one-dimensional  flow  initialization.  The  converged  coarse  grid 
solution  was  then  interpolated  onto  the  more  refined  grid,  thereby  providing  a  good 
initialization. 

The  convergence  history  for  the  full  scale  computation  is  shown  in  Figure  6.7. 

The  code  was  run  for  5000  iterations,  at  which  point  the  RMS  density  and  turbulent  kinetic 
energy  residuals  had  dropped  nearly  four  orders  of  magnitude.  (Note  that  the  total  number 
of  supersonic  points  (Mrei  >  0)  in  the  computational  domain  remained  nearly  unchanged  in 
early  iteration  due  to  the  coarse  grid  solution  initialization).  The  solution  required 
approximately  2.5  CPU  hours  on  a  Cray  Y-MP.  The  turbulence  intensity  was  arbitrarily 
set  to  3  %  ,  since  no  experimental  value  was  available,  and  the  artificial  dissipation  was 
scaled  by  the  reference  velocity  aooMrei  for  these  computations. 

In  Figure  6.8,  circumferentially  averaged  shroud  static  pressure  is  plotted  against 
time  averaged  static  tap  measurements  made  by  Krain  (1988).  The  agreement  is  very  good 
except  in  the  immediate  vicinity  of  the  impeller  inlet,  due  most  likely  to  the  cusped  leading 
edge,  which  is  relatively  thick  compared  to  the  blade  spacing  there  (refer  to  Figure  6.5). 


In  Figures  6.9  and  6.10,  predicted  meridional  velocities  (Wmejid  =  Vw£  +  W?) 


Iteration  # 


Figure  6.7  Convergence  history  for  Krain  impeller  How  computation. 
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Figure  6.8  Shroud  static  pressure  distribution.  Experimentally  obtained  time 
averaged  static  pressure  tap  measurements  (symbols),  circumferentially 
averaged  computed  values  (solid  line). 
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Figure  6.9  Meridional  velocity  profiles  at  L2F  data  acquisition  Plane  I. 
Experimental  measurements  (symbols),  computed  values  (solid  line). 
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are  compared  with  laser  measurements  at  planes  I  and  D,  which  correspond  to  0  %  and 
30  %  chord  as  measured  along  the  shroud  (refer  to  Figure  6.4).  Pitchwise  flow  turning  is 
significant  at  the  impeller  inlet,  since  the  blades  are  quite  thick  there  (elliptic  effect,  as  for  a 
turbine  blade).  The  resulting  blade-to-blade  pressure  gradient  gives  rise  to  larger 
meridional  velocities  near  the  blade  suction  side,  along  the  span.  This  is  captured  by  the 
simulation,  as  seen  in  Figures  6.9  and  6.10.  Rapid  axial  to  radial  turning  near  the  impeller 
inlet  gives  rise  to  a  positive  gradient  in  meridional  velocity  in  the  hub  to  shroud  direction  as 
seen  in  both  the  experimental  and  computational  results  in  Figure  6.10  (most  notable  if 
contrasted  with  inlet  profiles  in  Figure  6.9). 

In  Figure  6.11,  solution  contours  of  normalized  relative  helicity  are  plotted  at  a 
grid  slice  corresponding  to  Plane  HI  (48  %  chord).  In  this  figure,  high  positive  values  of 
normalized  helicity  (dark  shading)  correspond  to  a  stream  wise  vortex  (into  page),  whereas 
high  negative  values  correspond  to  negative  streamwise  vorticity.  The  secondary  motions 
arising  from  axial  to  radial  turning  of  the  blade  boundary  layers  (effect  A,  Figure  6.6a)  are 
clearly  evident.  Also,  the  tip  clearance  vortex  is  seen  to  have  migrated  well  into  the 
midpitch  region  near  the  shroud. 

Figure  6.12  shows  the  meridional  velocities  at  50  %  chord.  Both  experiment  and 
computation  clearly  show  a  region  of  significant  momentum  deficit  near  the  shroud.  This 
is  due  primarily  to  the  accumulation  of  low  momentum  boundary  layer  fluid  arising  from 
the  dominant  secondary  motions  discussed  previously.  Tip  clearance  fluid  is  also 
accumulated  in  this  "wake"  region.  It  was  found  that  the  level  of  agreement  between 
experiment  and  computation  in  this  region  was  only  achieved  when  the  number  of  grid 
points  in  the  tip  gap  was  increased  from  four  to  eight,  and  the  grid  was  hqi  clustered  in  the 
spanwise  direction  there  (refer  to  Figure  6.5c).  This  reinforces  that  both  axial  to  radial 
turning  induced  secondary  motions  and  tip  clearance  flow  play  an  important  role  in  the 
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Figure  6.12  Meridional  velocity  profiles  at  L2F  data  acquisition  Plane  m. 
Experimental  measurements  (symbols),  computed  values  (solid  line). 
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formation  of  this  low  momentum  fluid  region.  The  Navier-Stokes  solution  overpredicts  the 
extent  of  the  wake  region,  in  both  the  blade  to  blade  and  shroud  to  hub  directions.  In 
Figure  6.13,  meridional  velocity  contours  are  plotted.  Agreement  between  computation 
and  experiment  is  good,  although  it  appears  that  the  predicted  wake  flow  has  displaced  the 
high  momentum  inviscid  core,  or  "jet,"  fluid  towards  the  suction  surface  prematurely. 

Figure  6.14  shows  the  normalized  relative  helicity  contours  at  65  %  chord  (Plane 
TV).  The  jet-wake  structure  is  fully  formed  at  this  location.  Still  evident  are  th^  "ounter 
rotating  secondary  vortices.  The  tip  clearance  vortex  is  clearly  distinguishable,  and  appears 
to  have  entrained  some  of  the  fluid  from  the  pressure  side  half  of  the  passage,  "dragging" 
some  high  streamwise  vorticity  fluid  towards  the  suction  side.  Meridional  velocities 
compare  quite  well  at  Plane  IV,  as  shown  in  Figure  6.15.  Most  notable  is  that  the 
experimental  and  computed  wake  flow  has  extended  to  midspan.  The  blockage  that  this 
fluid  provides  is  seen  to  accelerate  the  higher  momentum  jet  fluid  towards  the  two  blade 
and  hub  surfaces.  Peak  values  of  meridional  velocity  occur  near  the  suction  side-hub 
comer.  These  features  are  evident  in  the  meridional  velocity  contours  provided  in  Figure 
6.16,  where  it  is  also  seen  that  the  location  of  low  meridional  velocity  near  the  shroud  has 
migrated  towards  the  suction  surface. 

The  results  at  65  %  chord,  presented  in  the  previous  paragraph,  illustrate  that  the 
wake  flow  has  fully  manifested  itself  at  this  location.  It  is  interesting  to  compare  a  highly 
referenced  qualitative  sketch  provided  by  Eckardt  (1976),  with  normalized  relative  helicity 
solution  contours.  In  Figure  6.17a,  Eckardt's  sketch,  which  serves  to  illustrate  the  variety 
of  phenomena  which  contribute  to  the  formation  of  the  jet-wake  flow,  has  been  adapted.  A 
solution  grid  slice  at  68  %  chord  was  chosen  for  comparison  in  Figure  6.17b,  since  the 
salient  features  of  the  wake  flow  correspond  well  with  Eckardt's  sketch  at  this  location. 
This  comparison  clearly  serves  to  illustrate  the  appropriateness  of  the  present  Navier- 
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Figure  6.15  Meridional  velocity  profiles  at  L2F  data  acquisition  Plane  IV. 
Experimental  measurements  (symbols),  computed  values  (solid  line). 
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Stokes  procedure  in  reconciling  these  flow  physics. 

At  80  %  chord  (Plane  V),  several  channel  vortices  are  clearly  distinguishable  as 
seen  in  Figure  6.18.  One  of  these,  of  positive  streamwise  vorticity,  has  migrated  from  the 
pressure  side  half  of  the  channel  (refer  to  Figures  6.11  and  6.14).  The  tip  clearance  vortex 
has  continued  to  migrate  closer  to  the  suction  surface.  The  meridional  velocity  profiles  in 
Figure  6.19  illustrate  that  only  the  near  hub  region  shows  inviscid  flow  features,  the  rest  of 
the  passage  flow  being  highly  distorted  by  the  variety  of  interacting  secondary  motions. 
This  is  captured  well  by  the  numerical  simulation,  as  seen  by  comparison  of  the 
experimental  and  computational  results  in  Figures  6.19  and  6.20. 

Near  the  impeller  discharge  plane  (Plane  VI,  98  %  chord),  the  entire  passage  flow 
is  encompassed  by  a  very  complex  shear  layer.  Figure  6.21  shows  normalized  relative 
helicity  contours  in  this  plane.  A  variety  of  nearly  indistinguishable  secondary  flow 
structures  are  present  Experimentally  measured  exit  velocities  are  rather  smooth,  though 
clearly  distorted  across  the  discharge  area.  This  is  illustrated  in  Figure  6.22. 

Computational  results  at  this  measurement  plane  are  poor,  due  to  the  "pinched"  trailing 
edge.  Meridional  velocity  profiles  for  the  discharge  plane  are  provided  in  Figure  6.23. 


Computed  and  measured  performance  parameters  are  presented  in  Table  6.1. 
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Figure  6. 19  Meridional  velocity  profiles  at  L2F  data  acquisition  Plane  V. 
Experimental  measurements  (symbols),  computed  values  (solid  line). 
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Figure  6.21  Computed  normalized  relative  helicity  contours  at  solution  grid 
slice  corresponding  to  L2F  data  acquisition  Plane  V. 


Figure  6.22  Meridional  velocity  profiles  at  L2F  data  acquisition  Plane  VI. 
Experimental  measurements  (symbols),  computed  values  (solid  line). 
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Table  6.1.  Comparison  of  Performance  Parameters  for 
Transonic  Centrifugal  Compressor  Stage 


In  order  to  obtain  some  qualitative  measure  of  the  importance  of  Reynolds  stress 
anisotropy  on  the  flowfield  structure  for  this  case,  and  to  examine  the  numerical  robustness 
of  the  hybrid  model  in  flows  where  significant  extra  strain  rates  are  present,  the  model  was 
used  to  compute  the  present  impeller  flowfield.  In  Figure  6.24,  a  convergence  history  for 
these  calculations  is  presented.  The  5000  iteration  k-e  solution  discussed  above  was  used 
to  initialize  the  ARSM  run  (compare  Figure  6.7).  The  jump  in  residuals  at  iteration  5000 
corresponds  to  a  solution  restart  with  the  new  model.  Specifically,  the  hybrid  model  was 
run  for  2880  iterations,  with  the  production  by  rotation  tensor,  Rjj  in  equation  [2.14],  set  to 
zero  (Q  =  0  in  equation  [3.60]).  This  solution  was  ran  until  the  RMS  density  residual 
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Figure  6.24  Convergence  history  for  the  three  computations  presented.  Features 
labelled  a  and  b  correspond  to  solution  restarts  using  alternate  turbulence  models. 
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recovered  to  a  four  decade  drop.  The  production  by  rotation  terms  were  then  implemented 
and  the  code  was  run  for  another  1480  iterations  where  again  the  log  of  the  density  residual 
had  dropped  to  -4.  All  three  solutions  are  compared  below.  It  is  noted  that  the  "peaks"  in 
density  and  turbulent  kinetic  energy  residual  which  develop  with  increased  iteration  count 
occur  at  100  iteration  intervals,  corresponding  to  solution  restart  frequency,  no!  the 
frequency  with  which  Reynolds  stresses  are  recomputed.  The  overall  convergence  rate  of 
the  ARSM  solutions  are  not  significantly  deteriorated,  and  it  is  hoped  that  the  mechanism 
which  causes  these  features  can  be  identified  and  fixed. 

In  Figure  6.25,  the  contour  of  y+match  =  200  for  the  hybrid  model  is  presented  for 
a  grid  slice  at  75  %  chord,  to  illustrate  the  proximity  of  the  hybrid  model  matching  location 
to  the  passage  surfaces. 

Figure  6.26  shows  converged  meridional  velocity  profiles  at  laser  Plane  V  for  the 
1)  k-e  (solid  line),  2)  hybrid,  Rjj  =  0  (long  dash),  3)  hybrid,  Ry  *  0  (short  dash)  solutions. 
The  meridional  velocity  profiles  are  seen  to  be  very  similar.  However,  the  near  wall  cross 
flow  velocities  and  wall  shear  stress  distributions  were  seen  to  vary  substantially  for  the 
various  models. 

Figures  6.27  shows  near  wall  cross  flow  velocity  profiles  near  the  exit  of  the 
impeller  (90  %  chord).  Figures  6.27a  and  6.27b  are  profiles  taken  at  midspan  on  the 
suction  and  pressure  surfaces  respectively,  whereas  Figure  6.27c  is  taken  at  mid-hub.  In 
these  plots,  the  square,  circular  and  triangular  symbols  correspond  respectively  to  the  k-e, 
hybrid,  Ry  =  0  and  hybrid,  Ry  *  0  solutions.  Qualitatively,  these  spanwise  flows  seem  to 
have  been  notably  influenced  by  turbulence  amplification  effects  on  the  developing  passage 
flow  boundary  layers,  as  seen  by  the  differences  in  the  cross  flow  velocity  profiles  in  these 
figures. 


1S9 


Figure  6.25  Converged  hybrid  model  solution  contours  of  on  a  grid  slice  at  75  % 
chord.  The  precise  location  of  y+match  =  200,  the  value  used  in  the  hybrid  solutions, 

is  indicated. 
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Figure  6.26  Meridional  velocity  profiles  at  L2F  data  acquisition  Plane  IV. 
Comparison  of  solutions  using  three  turbulence  models :  low  Reynolds  number 
k-e  model  (solid),  hybrid  model,  Ry  =  0  (long  dash),  hybrid  model,  Ry  *  0 

(short  dash). 
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Figure  6.27  Near  wall  relative  crossflow  velocity  profiles  at  90  %  chord, 
a)  midspan  on  suction  surface,  b)  midspan  on  pressure  surface,  c)  mid  passage 
on  hub.  Comparison  of  solutions  using  three  turbulence  models  :  low  Reynolds 
number  k-e  model  (square),  hybrid  model,  Rjj  =  0  (circle),  hybrid  model, 

Ry  *  0  (triangle). 
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In  Figure  6.28,  predicted  skin  friction  versus  chord  is  plotted  for  the  three 
solutions  at  midspan  on  the  suction  surface,  mid  span  on  the  pressure  surface,  and  at  mid¬ 
hub  locations.  In  this  figure,  the  solid  line  corresponds  to  the  k-E  model,  the  long  dash  to 
the  hybrid  model  without  rotation  terms,  and  the  short  dash  to  the  hybrid  model  with 
rotation  terms.  Choice  of  model  seems  to  have  little  effect  on  suction  surface  shear  stress. 
The  effect  of  extra  strain  rates  on  the  hybrid  model  do  seem  to  significantly  influence 
developing  boundary  layers  along  the  hub  and  pressure  surfaces. 

In  Figure  6.29,  normalized  relative  helicity  contours,  at  68  %  chord,  are  compared 
between  the  k-E  solution  and  the  hybrid  model  solution,  with  rotation  terms  included. 
Secondary  motions  are  qualitatively  very  similar  in  this  region,  where  the  jet-wake  flow  is 
fully  realized.  Quantitative  differences  are  certainly  noticeable.  These  results  do  not 
indicate  that  the  influence  of  curvature  and/or  rotation  on  turbulence  damping  in  the  shroud- 
suction  surface  region  significantly  influence  the  structure  of  the  wake. 

Table  6.2  compares  predicted  performance  parameters  between  the  three  modelling 
approaches. 


Figure  6.28  Predicted  skin  friction  coefficient  vs.  dimensionless  stream  wise 
distance  along  shroud.  Comparison  of  solutions  using  three  turbulence 
models  :  low  Reynolds  number  k-e  model  (solid),  hybrid  model,  Rjj  =  0  Gong 
dash),  hybrid  model,  Rjj  *  0  (short  dash). 
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Figure  6.29  Computed  normalized  relative  helicity  contours  on  a  solution  grid 
slice  at  68  %  chord.  Comparison  of  solutions  using  two  turbulence  models : 
low  Reynolds  number  k-e  model  (bottom),  hybrid  model,  Ry  *  0  (top). 


Table  6.2.  Comparison  of  Predicted  Performance  Parameters  for 
Transonic  Centrifugal  Compressor  Stage,  Using  Different  Turbulence 

Models 
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Mass  Averaged  Flowfield 

Computed 

Computed 

Computed 

Parameter3 

(k-E  model) 

(hybrid.Rj;  =  0) 

(hybrid,Rn  *  0) 

Pressure  Rise, 

PO,  diffuser  exit 

PO,  impeller  inlet 

3.80 

4.01 

3.92 

Slip  Factor, 

Vfl,  impeller  exit 

Uimpeller  exit 

0.85 

0.86 

0.86 

The  slip  factor,  which  arises  due  to  coriolis  "eddy"acceleration  even  in  the  absence 
of  shear  is  seen  to  be  negligibly  affected  by  turbulence  model  choice.  Predicted  stage  total 
pressure  rise  varies  by  as  much  as  5  %  with  different  turbulence  modelling. 

As  is  clear  from  the  foregoing  results,  the  flow  in  this  machine  is  very  complex. 
As  such,  it  is  difficult  in  this  case  to  make  quantitative  arguments  as  to  the  direct  influence 
of  the  turbulence  model  on  boundary  layer  structure.  Specifically,  it  would  be  too  bold  to 
conclude  that  local  extra  strains  in  the  algebraic  model  influence  local  boundary  layers 
direcdy  based  on  the  results  presented  here.  This  is  because  local  boundary  layer  physics, 
especially  near  the  exit,  are  strongly  dependent  on  the  history  of  the  entire  developing 
passage  flow. 

What  can  be  concluded  from  the  results  of  this  section  are  that  extra  strain  rates  are 
seen  to  influence  the  meridional  flow  physics  only  slightly,  through  local  turbulence 
structure  modification.  Secondary  motions,  wall  shear  stress  and  total  pressure  rise 
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performance  show  increas  ed  sensitivity  to  modelling  choice.  In  the  opinion  of  the  author, 
the  proposed  hybrid  k-e/ARS  model  is  a  sensible  engineering  approach  to  computing 
complex  flows  where  strong  extra  strains  are  present.  This  approach  is  computationally 
efficient,  reconciles  the  low  Reynolds  number  k-E  form  with  the  algebraic  stress  modelling, 
and  seems  to  be  numerically  stable.  However,  further  validation  work  is  still  needed, 
specifically  the  treatment  of  flows  with  isolated  extra  strain  rates. 
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CHAPTER  7 

CONCLUSIONS  AND  RECOMMENDATIONS 

A  new  three-dimensional  explicit  Navier-Stokes  procedure,  which  incorporates 
two-equation  turbulence  modelling  has  been  developed  and  presented  in  this  thesis. 

Several  analyses  and  techniques  which  have  enabled  convergent  and  accurate  simulation  of 
high  Reynolds  number  turbomachinery  flows  on  highly  stretched  grids  have  been 
developed  and  presented.  These  include  extensive  stability,  order  of  magnitude  and 
numerical  verification  studies  applied  to  the  discrete  system  of  seven  governing  equations, 
with  emphasis  on  turbulence  modelling.  A  compact  finite  difference  flux  evaluation 
procedure  was  presented.  Local  velocity  and  flux  Jacobian  scalings  of  artificial  dissipation 
operators  were  devised  and  applied.  A  hybrid  low  Reynolds  number  k-e/high  Reynolds 
number  algebraic  Reynolds  stress  model  was  developed.  The  procedure  has  been  applied 
to  a  variety  of  turbomachinery  flowfields,  across  a  wide  range  of  Mach  numbers.  The 
results  of  these  computations  were  presented  and  interpreted,  and  indicated  that  the  method 
provides  accurate  and  convergent  simulation  of  turbomachinery  flowfields. 

A  study  of  the  stability  of  the  coupled  system  of  seven  transport  equations, 
including  a  low  Reynolds  number  form  two-equation  turbulence  model,  was  presented. 
The  linear  analysis  procedures,  order  of  magnitude  analyses  and  subsequent  numerical 
verification  studies  presented,  provided  several  conclusions.  Most  notable  are  three  items 
which  directly  relate  to  the  turbulence  transport  modelling.  Specifically,  it  was  shown 
through  stability  and  order  of  magnitude  analysis,  and  numerical  verification  ,  that  the 
turbulence  model  source  terms  themselves  do  not  adversely  affect  the  stability  of  the 
explicit  numerical  procedure,  anywhere  in  the  flowfield,  except  in  early  stages  of  iteration. 
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This  is  a  result  contrary  to  that  generally  perceived.  Secondly,  it  was  shown  that 
turbulence  transport  models,  which  provide  large  values  of  effective  diffusion  in  regions 
where  grid  clustering  may  be  high,  give  rise  to  dominant  parabolic  stability  limitations, 
which  must  be  incorporated  in  the  numerical  procedure.  This  is  unlike  the  solely 
hyperbolic  constraints  often  sufficient  when  algebraic  eddy  viscosity  models  are  used  in 
time  marching  Navier-Stokes  procedures  (point  transition  specification  in  these  models 
preclude  large  values  of  eddy  viscosity  upstream  of  the  blade  row).  This  second  item  was 
clearly  demonstrated  as  a  mechanism  in  turbomachinery  blade  row  computations.  In  the 
opinion  of  the  author,  this  effect  must  manifest  itself  also  in  unstructured  grid 
implementations  where  high  local  clustering  appears  well  away  from  wall  damping  effects, 
as  for,  say,  shock  capturing.  Thirdly,  specific  near  wall  modelling  terms  which  appear  in 
the  stability  analysis  were  identified,  which  upon  inspection,  reveal  the  adverse  stability 
implications  of  incorporating  a  dissipation  rate  model  which  transports  a  total  dissipation 
rate  variable.  This  analysis  substantiates  claims  to  this  effect  which  appear  in  the  literature 
and  have  been  the  personal  experience  of  several  researchers,  including  the  author. 

Two  other  results,  related  to  the  numerical  implementation  of  the  turbulence  model, 
were  investigated  and  presented.  Specifically,  implicit  treatment  of  the  turbulence  source 
terms  did  not  improve  the  convergence  rate  of  the  explicit  scheme,  consistent  with  the  first 
conclusion  in  the  previous  paragraph.  Additionally,  terms  which  account  for  unstable, 
early  iteration  turbulence  model  behavior,  arising  from  inviscid  initialization  procedures, 
were  identified,  and  several  stabilization  procedures  that  the  author  has  found  to  be 
effective  were  presented. 

The  implementation  of  algebraic  and  full  Reynolds  stress  closures  in  complex 
aerodynamic  flow  simulations  should  increase  in  the  short  term.  This  will  be  due  to  the 
desire  to  move  to  improved  modelling,  increased  digital  computer  speed  and  memory,  and 


199 


the  onset  of  radically  different  designs  for  which  Navier-Stokes  methods  will  be  more 
critical  in  analysis.  The  turbulence  model  numerical  considerations  presented  in  Chapter  3, 
and  summarized  above,  should  be  useful  in  the  move  to  second  order  closures  for  two 
reasons.  Firstly,  the  modelled  FRS  and  dissipation  rate  equation  system  is  similar  to  the  k- 
e  system  in  that  they  are  both  convection-diffusion  transport  equations  with  large  non-linear 
source  terms.  Accordingly,  some  of  the  procedures  and  conclusions  presented  in  this 
thesis  will  be  applicable.  Secondly,  in  ARS  closures,  an  "anisotropic  eddy  viscosity"  can 
be  identified,  to  which  the  parabolic  stability  conclusions  recalled  above  are  directly 
applicable.  This  conclusion  was  demonstrated  analytically,  and  numerically  verified,  in 
Chapter  3. 

In  addition  to  turbulence  model  issues,  analyses  were  presented  which  focused  on 
the  influence  of  numerical  coupling  of  the  mean  flow  and  turbulence  transport  equations, 
system  rotation  source  terms  and  artificial  dissipation  on  the  numerical  stability  of  the 
Runge-Kutta  procedure.  Linear  analysis  and  numerical  verification  provided  that  there  is 
no  advantage  to  numerically  coupling  the  two-equation  model  system  to  the  mean  flow 
equation  system  in  regard  to  convergence,  accuracy  or  performance  of  the  scheme.  The 
appearance  of  rotation  source  terms  in  the  discrete  governing  equations  was  found  to  have 
negligible  influence  on  the  stability  of  the  scheme  for  rotation  numbers  and  grids  typical  of 
high  Reynolds  number  turbomachinery  rotor  computations.  It  was  found  useful  to  include 
artificial  dissipation  operators  in  the  specification  of  a  local  timestep  in  problems  where 
parabolic  stability  bounds  dominate  locally  (as  in  the  blade  row  computations  considered) 
to  avoid  ad  hoc  specification  of  an  operational  VonNeumann  number  which  varies  with 
level  of  artificial  dissipation  used. 

Some  of  the  limitations  of  the  foregoing  stability  results  were  also  presented.  In 
multigrid  procedures,  the  turbulence  model  and  rotation  source  terms  give  rise  to  stability 
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limitations  which  may  not  be  negligible  on  coarser  grids,  if  complete  residuals,  including 
source  terms,  are  evaluated  on  these  successively  coarser  grids.  Also,  the  flow  physics  are 
significantly  altered  by  system  rotation  .  Accordingly,  system  rotation  does  affect 
convergence  rates  since  stationary  and  rotating  flows,  through  the  same  geometry,  will  be 
very  different  (especially  in  centrifugal  machines).  The  related  conclusions  presented 
thereby  only  apply  to  the  specification  of  a  local  timestep.  Lastly,  the  analyses  presented 
are  local  and  linearized.  As  such,  the  conclusions  of  these  analyses  are  approximate,  and 
their  usefulness  had  to  be  verified,  like  any  engineering  approximation,  through  numerical 
application. 

Two  artificial  dissipation  issues  of  importance  in  high  Reynolds  number  flowfield 
computations  were  presented  and  discussed.  Two  scaling  functions  which  accommodate 
these  considerations  were  presented.  Specifically,  an  alternative  flux  Jacobian  eigenvalue 
scaling  was  introduced  and  extended  to  three  dimensions.  This  scaling  reduces  dissipation 
levels  in  the  directions  normal  to  grid  stretching  without  increasing  dissipation  levels  in  the 
same  direction  as  the  grid  stretching  (nominally,  the  wall  normal  direction).  It  was  found 
to  be  crucial  to  incorporate  this  anisotropic  eigenvalue  scaling  in  order  that  a  converged 
three-dimensional  rotor  flow  solution  be  obtained  on  a  highly  stretched  grid.  The  role  of 
artificial  dissipation  in  corrupting  solution  accuracy  in  near  wall  regions,  where  physical 
diffusion  is  important,  was  illustrated.  A  local  velocity  squared  scaling  which  tapers  levels 
of  artificial  dissipation  to  zero  near  solid  boundaries  was  introduced.  It  was  demonstrated 
through  numerical  experimentation  that  it  is  crucial  to  incorporate  this  scaling  in  order  that 
near  wall  physics  be  accurately  resolved. 


The  flux  evaluation  scheme  presented  eliminates  the  metric  ambiguity,  at  the 
interface  between  solid  and  periodic  boundaries,  associated  with  standard  finite  difference 
approaches.  This  scheme  allows  for  straightforward  inviscid  and  viscous  flux  evaluation, 
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at  complex  rotor  passage  boundaries,  including  periodic  cascade  and  tip  gap  boundaries, 
and  the  interface  between  rotating  and  stationary  sections  of  the  hub. 

A  hybrid  low  Reynolds  number  k-e/high  Reynolds  number  algebraic  Reynolds 
stress  model  was  developed.  This  model  reconciles  the  near  wall  damping  of  the  k-e 
model,  with  the  ability  of  the  ARS  model  to  provide  Reynolds  stress  anisotropies  which 
arise  due  to  extra  strain  rates.  The  formulation  of  this  model  and  its  numerical 
implementation  and  stability  topics  were  presented,  and  preliminary  model  validation  was 
provided  .  The  model  was  found  to  predict  near  wall  turbulence  behavior  correctly, 
including  normal  stress  anisotropy,  and  was  shown  to  be  stable  when  implemented 
numerically  with  a  local  blending  function.  This  model  is  computationally  efficient  in  full 
scale  three-dimensional  applications  because  the  influence  of  near  wall  pressure  strain 
terms,  which  are  very  complicated  in  three-dimensions,  are  implicitly  accommodated  in  the 
modelled  two-equation  near  wall  damping  terms.  In  the  opinion  of  the  author,  this  is  a 
sensible  engineering  approach  to  employing  a  second  order  closure  in  a  time-marching  full 
Navier-Stokes  procedure.  Specifically,  the  influence  of  curvature  and  rotation  on  Reynolds 
stress  anisotropy,  which  can  significantly  influence  mean  flow  structure  in  t-rbomachines 
is  incorporated,  but  wall  functions  are  not  used.  Therefore,  complex  near  wall  three- 
dimensional  boundary  layer  physics,  which  are  important  in  many  turbomachines  can  be 
explicitly  resolved  as  well.  Though  preliminary  validation  has  been  provided  and 
application  to  full  scale  turbomachinery  computations  demonstrated  in  this  thesis,  some 
further  study  may  be  useful.  Specifically,  the  application  of  the  model  to  configurations 
where  isolated  strains  give  rise  to  Reynolds  stress  and  mean  flow  redistributions  should  be 
further  investigated.  The  computation  of  a  simple  shear  layer  subject  to  spanwise  rotation 
is  under  way,  and  the  prediction  of  secondary  motion  of  the  second  kind  could  be 
demonstrated  for  a  square  duct  calculation. 
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Two-dimensional  results  were  presented.  Validation  of  the  turbulence  models  was 
provided,  and  the  results  of  one  supersonic  and  one  low  subsonic  compressor  cascade 
computations  were  presented.  Flowfield  predictions  were  found  to  be  acceptable  for  both 
cascades  except  in  the  aft  potion  of  the  subsonic  cascade  suction  surface,  due  to  mean  flow 
reversal  there.  Overall  cascade  performance  parameters  were  well  predicted  for  the 
supersonic  cascade  but  not  well  predicted  for  the  low  subsonic  cascade,  due  again  to 
flowfield  unsteadiness  as  well  as  turbulence  model  shortcomings. 

The  three-dimensional  procedure  was  used  to  predict  the  flow  in  a  compressor 
rotor.  The  viscous  incompressible  flow  through  a  curved  duct  was  predicted  accurately,  as 
verified  by  comparison  of  CFD  validation-quality  pressure  and  primary  and  secondary 
velocity  measurements  with  computed  values.  A  number  of  conclusions  apply  to  the  rotor 
flow  computation.  Specifically,  the  viscous  flow  through  a  rotor  was  predicted,  including 
inviscid  regions,  the  stagnation  region  near  the  leading  edge,  blade  boundary  layers  and 
wakes,  and  spanwise  mixing  region  downstream.  Flow  losses  and  loading  coefficient 
were  predicted  reasonably  well.  The  flow  in  the  tip  region  was  resolved  qualitatively,  with 
the  exception  of  nearwall  streamwise  and  radial  velocity  profiles.  The  effects  of  rotation 
and  endwall  secondary  flow  on  the  radial  velocity  profiles  in  the  near  wake  close  to  the  hub 
was  captured  reasonably  well.  The  ability  of  the  code  to  predict  spanwise  and  lateral 
mixing  downstream  of  the  rotor  was  demonstrated.  Specifically,  the  redistribution  of 
velocities  and  flow  angles  was  predicted  accurately.  Additionally,  the  spanwise 
redistribution  of  losses  between  the  trailing  edge  and  the  far  wake  region  was  captured 
well. 


The  flowfield  in  a  transonic  centrifugal  compressor  stage  (impeller  +  diffuser)  was 
computed.  The  numerical  and  modelling  methodologies  provided  in  this  thesis  give  rise  to 
a  reasonably  accurate  simulation  of  this  complex  flowfield.  Specifically,  stage  pressure 
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rise  distribution  was  captured  accurately.  Also,  the  secondary  motions  arising  from 
coriolis  and  centrifugal  effects  and  tip  clearance  flow  were  resolved  with  reasonable 
accuracy.  Normalized  relative  helicity  was  helpful  in  qualifying  and  tracking  some  of  these 
motions.  The  accumulation  of  low  momentum  fluid  near  the  shroud,  in  the  aft  portion  of 
the  impeller  was  captured.  This  fluid  includes  that  convected  from  blade  boundary  layers 
by  curvature  induced  secondary  motion  and  tip  clearance  flow,  and  accumulates  in  the  mid¬ 
passage  region.  Near  the  impeller  exit,  this  fluid  extends  well  into  the  impeller  passage 
towards  the  hub,  causing  blockage  (and  hence  acceleration)  of  the  high  momentum  fluid 
near  the  suction  surface-hub  comer.  These  physics  are  evident  in  experimental  L2F 
meridional  velocity  measurements  and  in  the  numerical  simulation. 

The  hybrid  k-e/ARS  model  was  applied  to  this  flowfield  in  order  to  obtain 
qualitative  measure  of  the  influence  of  Reynolds  stress  anisotropy  on  the  developing 
flowfield  structure  in  the  centrifugal  impeller,  and  to  examine  the  numerical  robustness  of 
the  model  in  flows  where  significant  extra  strain  rates  are  present.  This  is  the  highest  level 
of  turbulence  modelling  ever  attempted,  to  the  knowledge  of  the  author,  in  this  type  of 
machine.  It  was  found  that  extra  strain  rates  influence  meridional  flow  physics  and 
impeller  slip  factor  only  slightly,  through  local  turbulence  structure  modification.  Cross 
flows,  wall  shear  stress  and  total  pressure  rise  performance  showed  slightly  increased 
sensitivity  to  modelling  choice.  The  model  was  found  to  be  stable  for  this  computation, 
though  solution  restart  anomalies  remain  to  be  resolved. 

In  light  of  the  relatively  small  difference  in  predicted  meridional  physics  arising 
due  to  modelling  choice,  when  compared  to  the  differences  in  computed  and  measured 
meridional  flow,  geometric  modelling  considerations  are  a  priority,  in  order  to  provide 
improved  prediction  of  these  types  of  flows.  Specifically,  inclusion  of  the  actual  spinner 
geometry,  though  difficult  from  a  grid  topology  standpoint,  should  provide  more  accurate 
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inlet  hub  boundary  layer  resolution  (and  thereby  the  influence  of  the  developing  secondary 
motions  which  arise  due  to  the  turning  and  rotation  of  this  shear  layer.)  More  important, 
though,  is  the  accurate  resolution  of  the  leading  edge,  trailing  edge  and  tip  gap  regions. 
The  "pinched"  H-grid  approach  at  these  boundaries  is  most  likely  the  dominant  source  of 
error  in  these  computations,  since  blades  are  thick  compared  to  passage  dimensions,  and 
the  tip  clearance  flow  is  so  important  in  centrifugal  impellers.  Lastly,  as  mentioned  above, 
further  validation  is  required  for  the  turbulence  models  used,  before  stronger  conclusions 
can  be  drawn  in  regard  to  application  of  the  model  to  a  flow  as  complex  as  that  in  a 
centrifugal  compressor. 
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APPENDIX 

NONCONSERVATIVE  JACOBIAN  MATRICES  AND 
FOURIER  MATRIX  FOR  HYPERBOLIC  STABILITY 


Nonconservative  Jacobian  Matrices 


Fourier  Matrix  for  Hyperbolic  Stability 
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